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ABSTRACT 


We  present  a  distributed,  physically— based  model  of  runoff  generation  in  a 
catchment,  for  operational  use  in  flood  forecasting.  The  model  accounts  for  both  the 
infiltration— excess  and  saturation-excess  mechanisms  of  runoff  production,  and  for 
lateral  subsurface  flows.  Tie  effect  of  local  terrain  slope  on  subsurface  flows  and  the 
development  of  areas  of  saturated  soil  is  accounted  for.  The  model  uses  spatial 
discretization  into  rectangular  elements  which  correspond  to  the  grid  of  a  Digital 
Elevation  Map.  Each  basin  element  consists  of  a  soil  column  in  which  hydraulic 
conductivity  decreases  with  depth,  in  the  form  of  an  exponential  function.  Spatial 
discretization  allows  for  distributed  terrain  slope,  soil  parameters,  moisture 
conditions,  and  rainfall  inputs.  Time  discretization  allows  for  consideration  of 
time-variable  rainfall  rates. 

The  model  uses  the  kinematic  approximation  of  infiltration  and  subsurface 
flows,  which  are  assumed  to  occur  only  within  the  porous  soil  matrix.  The 
kinematic  model  of  infiltration  is  used  to  show  how  decreasing  conductivity  with 
depth  may  result  in  the  development  of  a  zone  of  perched  saturation  during  a 
rainstorm,  and  that  the  flow  in  the  perched  saturated  zone  is  diverted  laterally  ,f 
the  terrain  is  inclined.  A  simplified  computational  procedure  is  introduced  that 
allows  flow  transfer  among  elements.  Inter-element  flow  transfers  are  used  to 
predict  the  position  of  the  permanent  water-table  in  each  element  at  the  time  of 
initiation  of  the  rainfall  event  and  the  extent  of  the  area  of  saturated  soils. 

Moisture  transfer  between  elements  during  the  storm  is  also  considered.  The  model 
is  extended  to  consider  anisotropic  soils,  and  we  show  that  higher  lateral  than 
vertical  conductivity  results  in  increased  lateral  diversion  of  flow. 

The  model  was  applied  to  the  Sieve  catchment  in  Italy  and  used  to  reproduce 
hydrographs  for  12  recorded  rainstorms.  Given  that  pre-storm  basefiow  was  not 
available  for  any  of  the  12  storms,  three  different  water-table  positions  were 
considered  for  each  storm,  low  average  and  high,  corresponding  to  steady-state 
equilibrium  with  the  basefiow  values  that  have  a  90%  and  10%  probability  of  being 
exceeded  in  the  month  in  question.  The  observed  hydrographs  for  the  various 
storms  were,  generally,  in  the  area  comprised  between  the  “dry  ”  and  Vet  ” 
predicted  hydrographs,  or  did  not  fall  far  from  this  area.  Given  the  scarcity  of  data 
available  for  the  Sieve,  we  consider  model  predictions  to  be  quite  encouraging. 
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an  average  month  of  November.  The  water-table  is  in 
steady-state  with  the  basin  discharge  rate,  Qb  =  4.0  m3/s 

that  has  an  exceedence  probability  of  50%  in  November.  140 
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4.12  Depth  below  the  surface  (in  meters)  of  the  water-table  in 
a  wet  month  of  November.  The  water-table  is  in  steady- 
state  with  the  basin  discharge  rate,  =  10.0  m3/s  that 

has  an  exceedence  probability  of  10%  in  November.  141 

4.13  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1968. 
Groundwater  was  initialized  in  steady-state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  February,  Rj  =  0.0129 

mm/hr.  142 

4.14  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1968. 
Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  February, 

Rj  =:  0.0386  mm/hr.  143 

4.15  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  February,  1968. 

Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  February,  Rj  =  0.0728 

mm/hr.  144 

4.16  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1968. 

Groundwater  was  initialized  in  steady-state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  December, 

Rj  =  0.0087  mm/hr.  145 

4.17  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1968. 

Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  December, 

Rj  =  0.0300  mm/hr.  146 
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4.18  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1968. 
Groundwater  was  initialized  in  steady- state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  December, 

Rj  =0.0664  mm/hr.  147 

4.19  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1969. 
Groundwater  was  initialized  in  steady-state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  January,  Rj  =  0.0129 

mm/hr.  148 

4.20  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  January,  1969. 

Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  January, 

Rj  =  0.0321  mm/hr.  149 

4.21  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  January,  1969. 

Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  January, 

Rj  =  0.0171  mm/hr.  150 

4.22  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1975. 

Groundwater  was  initialized  in  steady-state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  December, 

Rj  =  0.0087  mm/hr.  151 

4.23  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1975. 

Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
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probability  of  50%)  in  the  month  of  December, 

Rj  =  0.0300  mm/hr.  152 

4.24  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1975. 
Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  December, 

Rj  =  0.0664  mm/hr.  153 

4.25  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1976. 
Groundwater  was  initialized  in  steady- state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  December, 

Ri  =  0.0087  mm/hr.  154 

4.26  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1976. 
Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  December, 

Rj  =  0.0300  mm/hr.  155 

4.27  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1976. 
Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  December, 

Rj  ss  0.0664  mm/hr.  156 

4.28  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1977. 
Groundwater  was  initialized  in  steady-state  with  the 
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dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  February,  Rj  =  0.0129 

mm/hr.  157 

4.29  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1977. 
Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  February,  Rj  =  0.0386 

mm/hr.  158 

4.30  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  cf  February,  1977. 
Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  February,  Rj  =  0.0043 

mm/hr.  159 

4.31  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1979. 
Groundwater  was  initialized  in  steady- state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  January,  Rj  =  0.0129 

mm/hr.  160 

4.32  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1979. 
Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  January,  Rj  =  0.0321 

mm/hr.  161 

4.33  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1979. 
Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
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probability  of  10%)  in  the  month  of  January,  Rj  =  0.0643 

mm/hr.  162 

4.34  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1981. 
Groundwater  was  initialized  in  steady- state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  December, 

Rj  =  0.0087  mm/hr.  163 

4.35  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1981. 
Groundwater  was  initialized  in  steady- state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  December, 

Rj  =  0.0300  mm/hr.  164 

4.36  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1981. 
Groundwater  was  initialized  in  steady- state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  December, 

Rj  =  0.0664  mm/hr.  165 

4.37  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  November,  1982. 
Groundwater  was  initialized  in  steady-state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  November, 

Rj  =  0.0043  mm/hr.  166 

4.38  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  November,  1982. 
Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  November, 

Rj  =  0.0171  mm/hr.  167 
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4.39  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  November,  1982. 
Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  November, 

R*  =  0.0428  mm/hr.  168 

4.40  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1983. 
Groundwater  was  initialized  in  steady-state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  February, 

Rj  =  0.0129  mm/hr.  169 

4.41  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  February,  1983. 

Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  February, 

Rj  =  0386  mm/hr.  170 

4.42  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  February,  1983. 

Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  February, 

Rj  =  0.0043  mm/hr.  171 

4.43  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  January,  1985. 

Groundwater  was  initialized  in  steady-state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  January,  Rj  =  0.0129 

mm/hr.  172 

4.44  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  January,  1985. 

Groundwater  was  initialized  in  steady-state  with  the 
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average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  January, 

Rj  =  0.0321  mm/hr.  173 

4.45  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1985. 
Groundwater  was  initialized  in  steady-state  with  the 
wet  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  10%)  in  the  month  of  January, 

Rj  =  0.0643  mm/hr.  174 

4.46  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  November,  1987. 
Groundwater  was  initialized  in  steady- state  with  the 
dry  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  90%)  in  the  month  of  November, 

Rj  =  0.0043  mm/hr.  175 

4.47  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  November,  1987. 
Groundwater  was  initialized  in  steady-state  with  the 
average  recharge  rate  (i.  e.  that  has  an  exceedence 
probability  of  50%)  in  the  month  of  November 

Rj  =  0.0171  mm/hr.  176 

4.48  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  November,  1987. 
Groundwater  was  initialized  in  steady-state  with  the 
probability  of  10%)  in  the  month  of  November, 

Rj  =  0.0428  mm/hr.  177 

4.49  Deviation  of  predicted  from  observed  peak-flow  versus 
cumulative  precipitation  in  the  30  days  preceding  the 

storm.  178 

A.l  Components  of  flow  in  the  z  and  h,  and  p  and  n 
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Chapter  1 


INTRODUCTION 


1 A.  STATEMENT  OF  THE  PROBLEM 


The  history  of  physically-based  modeling  of  runoff  generation,  no 
more  than  30  years  old,  has  been  characterized  by  an  enormously  wide 
range  of  approaches  and  levels  of  complexity.  While  part  of  the 
explanation  for  this  variety  lies  in  the  multiplicity  and  complexity  of  the 
physical  processes  themselves,  much  of  the  history  of  these  models  has 
been  shaped  by  the  gradual  development  of  computational  capabilities  and 
data  sources,  and  the  continuous  attempt  to  make  the  best  possible  use  of 
the  new  resources. 

One  of  the  major  difficulties  in  understanding  and  quantifying 
runoff  generation  in  river  basins  stems  from  the  presence  of  spatial 
variability  in  topography,  geology,  soil  type,  vegetation,  etc.,  and  in  water 
fluxes,  such  as  rainfall,  infiltration  and  evapotranspiration.  The 
accelerated  growth  of  remotely-sensed  distributed  data,  including  satellite 
readings  of  terrain  topography,  soil  type,  soil  use,  vegetational  coverage, 
geology,  soil  moisture  status,  and  radar  estimates  of  rainfall  rates,  brings 
enormous  promise  to  hydrologic  science  and  runoff  forecasting  in 
particular. 
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Such  distributed  data  are  rapidly  becoming  widespread  and  easily 
available  over  vast  regions  of  the  World.  In  particular,  Digital  Elevation 
Maps  (DEMs)  are  already  available  for  most  of  the  U.S.  territory,  and  other 
regions  of  the  World.  Weather  radar  provides  high  temporal  and  spatial 
resolution  estimates  of  rainfall-rates  and  benefits  rainfall  prediction 
enormously  through  detection  of  storm-system  development.  Near-full 
coverage  of  the  U.S.  and  Western  Europe  by  weather  radar  networks  is 
under  implementation.  Therefore,  there  is  a  challenge  to  develop  new 
model  formulations  which  are  capable  of  making  best  use  of  the 
distributed  data  available. 

The  availability  of  Digital  Elevation  Maps  (DEMs)  have  led  to  the 
generation  of  algorithms  for  estimation  of  topographic  variables  such  as 
terrain  slopes  and  exposure,  and  to  more  involved  ones  aiming  at  the 
derivation  of  river-channel  networks  from  a  number  of  assumptions  (e.  g. 
Tarboton,  1989).  Through  terrain  slopes,  the  runoff  model  is  able  to 
incorporate  the  effects  of  topography,  which  have  long  been  recognized  as 
important  in  runoff  generation.  Computer-derived  channel  networks  have 
revealed  a  generally  good  correspondence  with  map-digitized  networks, 
and  represent  a  convenient  alternative  to  them  for  use  in  distributed 
rainfall-runoff  models. 

Weather  radar  offers  unprecedented  perspectives  to  flood 
forecasting,  in  two  major  ways.  First,  it  provides  rainfall  rate  estimates 
with  high  spatial  and  temporal  resolutions.  While  such  estimates  may  not 
have  the  precision  of  raingage  point  measurements,  they  provide  the 
information  on  spatial  distribution  of  rainfall  that  is  crucial  to  accurate 
streamflow  forecasting.  Radar  makes  millions  of  measurements  per 
minute  with  a  resolution  in  the  order  of  one  kilometer,  and  its  readings 
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are  immediately  available.  Second,  rainfall  prediction  benefits 
enormously  from  radar  data.  Readings  at  successive  time  steps  reveal 
storm-system  development  and  movement.  The  more  advanced  Doppler 
technology  radar  offers  additional  advantages  as  a  result  of  its  ability  to 
detect  direction  and  speed  of  object  movement  in  a  single  reading.  A  color- 
code  map  of  air-mass  directions  and  speeds  can  be  produced  on  a 
computer  screen  in  real  time,  to  be  promptly  interpreted.  The  weather 
forecaster  can  thus  detect  rainstorm  development  and  movement,  and 
identify  the  particular  color  patterns  that  are  signatures  of  tornado 
development. 

Although  physical  arguments  favour  an  understanding  of  hillslope 
hydrology  on  the  basis  of  flow  strips,  remotely  sensed  data  are  generally 
available  for  grid  squares.  Both  DEMs  and  radar  estimates  of  rainfall  are 
made  available  for  grid  squares.  Growth  of  Geographic  Information 
Systems  data  bases  also  promise  to  provide  soil-type  maps  in  digital  form, 
from  which  parameter  values,  such  of  saturated  hydraulic  conductivity, 
may  be  derived.  Therefore,  there  seems  to  be  an  advantage  to  developing 
rainfall-runoff  models  structured  to  make  use  of  distributed  data  provided 
for  grid  squares.  From  a  computational  point  of  view,  spatial 
discretization  into  grid  squares  rather  than  flow  strips  or  irregular 
elements  facilitates  element  reference  through  coordinate  values, 
algorithm  structure,  and  data  storage  in  a  simple  matrix  form.  Also, 
algorithms  written  for  square  elements  are  easily  transferrable  for 
application  in  different  catchments. 

We  present  a  model  framework  that  incorporates  a  variety  of 
distributed  data  characterizing  terrain  morphology,  soil  parameters,  and 
rainfall  inputs,  in  the  form  of  grid  squares.  Terrain  morphology  is 
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obtained  from  DEMs,  soil  parameters  are  obtained  from  digitized  maps  of 
soil  type,  and  rainfall  rate  may  be  provided  by  radar  readings  or  by  a 
network  of  telemetering  raingages.  The  model  attempts  to  describe  the 
various  physical  processes  of  runoff  generation  on  the  hillslope  that  are 
currently  accepted,  namely  saturation-excess  runoff,  infiltration-excess 
runoff,  interflow,  and  return  flow. 


LB  REVIEW  OF  LITERATURE 


A  few  models  have  been  presented  which  seek  full  numerical 
solutions  to  the  governing  equations  of  subsurface  and  overland  flow. 
While  such  models  have  contributed  to  our  understanding  of  runoff 
production  processes  at  the  hillslope  scale,  their  data  and  computational 
requirements  precludes  their  application  at  the  catchment  scale.  One  of 
the  most  complete  of  such  models  is  that  of  Freeze  (1971),  which  simulates 
transient  subsurface  flow  in  an  integrated  saturated-unsaturated  flow 
system.  The  governing  equations  are  (1)  the  two-dimensional  equation  for 
Darcian  flow  and  (2)  the  Saint-Venant  equation  for  overland  flow.  The 
Darcian  flow  equation  uses  values  of  conductivity  and  hydraulic  potential 
measured  in  the  field.  Soil  conditions  of  heterogeneity  and  anisotropy  may 
be  considered  by  the  model.  Both  infiltration-excess  (Hortonian)  and 
saturation-excess  runoff  generation  are  predicted  by  this  model.  The 
model  admits  any  generalized  configuration  of  all  pertinent  boundary 
conditions,  namely  flow  divides,  position  of  impermeable  bedrock,  stream 
outflow  and  rainfall  input.  This  model  has  been  applied  to  a  hillslope  flow 
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strip  in  the  Reynolds  Creek  Watershed,  Idaho  (Stephenson,  G.R.  and 
Freeze,  R.A.,  1974).  Model  predictions  of  water-table  position,  hydraulic 
head  patterns,  and  generated  runoff  compared  well  with  field 
measurements  on  a  hillslope  flow  strip,  giving  support  to  the  model's 
construction.  However,  and  in  spite  of  the  small  scale  of  application  and 
the  large  data  collection  efforts,  the  model's  intense  requirements  for  field 
data  restricted  the  use  of  some  of  its  capabilities  such  as  consideration  of 
soil  anisotropy,  hysteresis,  and  of  streamflow  boundary  conditions. 

Model  application  at  the  catchment  scale  requires  introduction  of 
approximations  to  reduce  data  and  computational  requirements.  A 
generalized  approximation,  supported  by  field  evidence,  considers 
unsaturated  subsurface  flows  to  be  vertical  and  saturated  flows  to  be 
parallel  to  the  terrain  surface.  Saturated  flows  are  usually  treated  as  a 
single  layer,  which  reduces  the  problem  to  a  series  of  coupled  one¬ 
dimensional  equations  (of  similar  form  to  the  Saint-Venant  equations  for 
overland  flow).  One  model  that  considers  such  simplifications  is  the 
"Syst&me  Hydrologique  Europgen  (SHE)"  (Beven,  K.  et  al.,  1980;  Abbot, 
M.B.  et  al.,  1986a, b;  Bathurst,  1986a, b),  developed  through  the  combined 
efforts  of  the  Institute  of  Hydrology  (U.K),  the  SOGREAH,  and  the  Danish 
Hydraulic  Institute.  Some  models  that  include  these  simplifications  use 
the  kinematic  approximation  to  model  downslope  saturated  flow. 
Examples  are  the  model  of  Beven  (1981,  1982)  and  the  model  of  Hurley  and 
Pantelis  (1985). 

Some  runoff  models  emphasize  saturation-excess  over  infiltration- 
excess  runoff  generation.  In  predicting  saturation-excess  runoff,  the 
position  of  the  water-table  is  of  concern.  Several  researchers  have 
demonstrated  that  water-table  position  is  strongly  associated  with  terrain 


27 


morphology.  Aside  from  any  effect  due  to  random  heterogeneity  of  soil 
properties,  there  is  a  general  tendency  for  lateral  increase  in  moisture 
content  in  the  downslope  direction  that  results  from  the  increasing  size  of 
the  upstream  hillslope  area  contributing  subsurface  flows,  especially  in 
areas  of  convergent  topography.  If  the  upstream  drained  area  is 
sufficiently  large,  the  lower  parts  of  the  hillslopes  may  be  completely 
saturated,  or  water-logged.  The  significance  of  these  saturated  zones  is 
well  appreciated  in  disciplines  that  include  engineering,  soil  conservation, 
forestry,  agriculture,  and  hydrology  because  of  their  enormous  effect  on 
soil  strength,  erosion,  secondary  salinization,  and  storm  runoff.  Field 
studies  that  have  validated  the  importance  of  topography  in  determining 
the  position  and  extent  of  saturated  areas  include  the  pioneering  work  by 
Dunne  and  Black  (1970a,b),  and  the  studies  by  Anderson  and  Burt  (1977, 
1978)  and  Beven  (1978). 

Given  the  association  between  topography  and  water-table  position, 
many  models  that  emphasize  saturation-excess  runoff  generation  account 
for  topographic  variables  explicitly.  Dunne  and  Black  (1970)  were  the  first 
to  demonstrate  how  areas  of  saturated  soil  occur  most  readily  in 
convergent,  concave  areas,  and  that  these  areas  expand  and  contract 
during  rain  storms.  Similar  patterns  have  been  demonstrated  for 
subsurface  flows  by  Anderson  and  Burt  (1978),  although  Anderson  and 
Kneale  (1982)  have  shown  that  in  areas  of  low  slope  there  is  not  such  a 
close  relationship  between  saturation  and  surface  morphology.  Engman 
and  Rogowski  (1974)  combined  infiltration  equations  with  mapped  soil 
properties  to  forecast  saturated  areas.  Troendle  (Troendle,  1985;  Hewlett 
and  Troendle,  1975)  has  developed  a  series  of  related  models  which  use  a 
two-dimensional  numerical  solution  for  saturated  and  unsaturated  water 
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movement  along  a  flow  strip.  Where  surface  elements  near  the  slope  base 
become  saturated,  the  saturated  area  is  added  to  the  effective  channel  area 
receiving  direct  precipitation  and  the  remainder  of  the  hillside  upslope  is 
redivided  into  elements,  with  their  highest  density  close  to  the  edge  of  the 
saturated  area.  The  model  is  designed  to  give  its  most  accurate  forecasts 
for  the  margins  of  the  saturated  area,  which  then  generates  overland  flow. 
The  model  assumes  that  no  Hortonian  overland  flow  occurs.  The  model 
gave  good  agreement  for  a  0.38  Km2  mountain  catchment  in  West 
Virginia.  Ap^hcation  to  the  0.24  Km2  Whitefall  watershed  near  Athens, 
Georgia  was  less  successful. 

The  generalized  unavailability  of  field  data  with  the  detail  required 
for  physically-based  modeling  has  generated  strong  motivation  for 
development  of  conceptual  models.  Two  conceptual  models  have  been 
proposed  which  consider  catchment  morphology  to  predict  areas  of 
saturated  soils.  A  central  assumption  of  these  models  is  the  hypothetical 
condition  of  equilibrium  of  the  groundwater  storage  with  recharge.  One 
model  is  TOPMODEL  of  Beven  et  al.  (Beven  and  Kirkby,  1979;  Beven  and 
Wood,  1983;  Beven  et  al.,  1984;  Horberger  et  al.,  1985),  which  considers  an 
exponential  decrease  of  soil  hydraulic  conductivity  with  depth,  from  which 
results  an  exponential  form  for  the  relation  between  transmissivity  and 
water  table  depth.  The  model  has  been  extended  by  Beven  (1986a, b)  to  allow 
distributed  soil  parameters  and  to  include  prediction  of  infiltration-excess 
runoff.  TOPMODEL  was  able  to  give  moderately  good  forecasts  over  both 
storm  and  inter-storm  periods  in  most  applications.  The  other  model,  by 
O'Loughlin  and  Moore  (O'Loughlin,  1981,  1986;  Moore  et  al.,  1986)  can 
consider  arbitrary  variations  of  hydraulic  conductivity  with  depth.  Each  of 
these  two  models  involves  a  topographic  index  for  a  large  number  of  points 
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in  the  catchment,  and  the  two  indices  are  basically  very  similar  to  each 
other. 

The  topographic  index  introduced  in  TOPMODEL,  ln(a/tan(j3)),  where 
a  is  the  area  drained  per  unit  contour  length  and  tan((3)  is  the  slope  of  the 
ground  surface  at  the  location,  has  been  found  by  Beven  and  Kirkby  to 
compare  favorably  with  observed  patterns  of  surface  saturation.  Beven 
and  Kirkby  developed  a  hydrological  forecasting  model  combining  a 
dynamic  contributing  area  model  based  on  ln(a/tan((3)),  the  channel 
network  topology,  and  a  simple  lumped  parameter  basin  routing  model. 
The  model  parameters  were  physically  based  in  the  sense  that  they  may  be 
determined  directly  by  measurement.  Using  only  estimated  and 
measured  parameter  values  the  model  made  satisfactory  predictions  of 
catchment  response. 

O'Loughlin  (1981)  derived  criteria  for  the  existence  of  saturated 
areas  on  draining  hillslopes  in  natural  catchments,  based  on  the 
correlation  between  seepage  face  and  contributing  area.  These  criteria 
were  expressed  in  terms  of  a  topographic  index  h(x),  which  turns  out  to  be 
almost  identical  to  ln(a/tan((5)).  Based  on  the  analysis  of  hillslopes  of 
different  shapes  (parallel,  converging,  and  diverging),  O'Loughlin  showed 
that  the  contributing  areas  are  largest  in  the  convergent  zones  and 
smallest  in  the  divergent  zones.  Also,  the  uphill  boundary  of  the 
contributing  area  on  a  convergent  hillslope  was  shown  to  be  much  more 
stable  than  that  on  a  plane  or  divergent  hillslope.  This  means  that  in  drier 
conditions,  the  seepage  face  on  a  plane  hillslope  may  shrink  back  to  the 
stream  edge  while  on  a  convergent  slope  the  contributing  area  tends  to 
persist.  Earlier,  Kirkby  and  Chorley  (1967)  had  shown  theoretically  that 
the  largest  contributions  to  runoff  came  from  hillslopes  which  are  concave 
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in  section  and  convergent  in  plan.  These  predictions  have  been 
qualitatively  confirmed  by  the  observations  of  Dunne  (1978)  and  others  in 
the  field. 


1.C  MODEL  CONCEPTUALIZATION 


We  present  a  model  of  runoff  production  for  real-time  streamflow 
forecasting  which  considers  spatially  distributed  soil  parameters,  terrain 
slope,  and  rainfall  inputs.  The  model  is  physically  based  and  accounts  for 
both  infiltration-excess  (Horton  type)  and  saturation-excess  (Dunne  type) 
mechanisms  of  runoff  production.  The  kinematic  approximation  is  used 
to  model  subsurface  flows,  and  to  predict  infiltration  and  lateral 
groundwater  flows.  The  effect  of  local  terrain  slope  upon  lateral 
subsurface  flows  is  accounted  for.  Local  slope  is  obtained  from  a  Digital 
Elevation  Map  (DEM)  of  the  catchment,  which  is  provided  in  the  form  of 
elevation  values  taken  at  regularly-spaced  points  of  a  rectangular  grid. 
All  runoff  generated  is  routed  over  the  hillslope  to  a  channel  and  along  the 
channel  network  to  the  outlet,  and  contributes  to  the  predicted  hydrograph. 
Routing  considers  two  velocity  parameters  that  are  spatially  uniform  and 
constant  in  time;  a  velocity  of  flow  travel  over  hillslopes,  and  a  velocity  of 
travel  in  the  channel  network. 

To  account  for  spatial  variability  of  parameters  and  variables,  the 
basin  is  subdivided  into  identical  rectangular  elements,  or  pixels 
(contraction  of  the  term  picture  element).  The  use  of  rectangular  elements 
has  considerable  advantages  regarding  computational  efficiency  and  data 
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storage.  The  elements  are  taken  to  correspond  to  the  DEM  grid.  The  grid 
size  of  the  DEM  is  generally  comparable  to  the  resolution  of  digitized  soil 
maps  and  of  radar  rainfall  measurements,  so  that  distributed  information 
is  preserved  by  this  element  size.  Terrain  slope  for  each  element  is 
obtained  from  the  DEM.  The  network  of  stream  channels  may  be  digitized 
from  a  terrain  map  or  may  be  approximated  by  a  network  estimated 
automatically  from  the  DEM,  through  well  established  algorithms 
(Tarboton,  1989).  Automatically-generated  river  networks  have  in  general 
compared  favorably  with  mapped  blue  lines. 

Figure  1.1  represents  an  example  of  the  spatial  discretization  of  a 
river  basin  corresponding  to  a  DEM  in  a  rectangular  (in  this  example  it  is 
square)  grid. 


LD  OUTLINE 


This  work  is  organized  into  four  main  chapters,  followed  by  a 
concluding  chapter  and  an  appendix  including  the  computer  codes  used. 

Chapter  2  describes  the  model  of  infiltration,  which  is  based  in  the 
kinematic  approximation.  The  chapter  is  divided  into  three  sections. 
Section  2.A  provides  a  review  of  the  kinematic  model  of  infiltration  in  a 
homogeneous,  isotropic  soil,  following  the  presentation  by  Charbeneau 
(1984).  Section  2.B  presents  an  analysis  of  kinematic  infiltration  in  a 
vertically  heterogeneous  sloped  soil,  where  hydraulic  conductivity 
decreases  exponentially  with  depth.  Section  2.C  introduces  a  conceptual 
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approximation  to  the  a  complex  moisture  profile  which  allows 
computational  simplicity  in  the  treatment  of  variable  rainfall  rates. 

Chapter  3  describes  the  conceptualization  of  the  rainfall-runoff 
model  at  the  basin  scale.  The  chapter  is  divided  into  six  sections.  Section 
3.A  discusses  spatial  discretization  and  describes  computations  in  one 
element.  Section  3.B  discusses  the  computational  procedure  that 
accomplishes  element  coupling  through  lateral  flows.  Section  3.D  treats 
model  initialization  with  an  estimated  water-table  and  unsaturated 
moisture  profile  above  the  water-table.  Section  3.E  describes  the  flow 
routing  computation.  Finally,  Section  3.F  investigates  model  sensitivity  to 
soil  parameters,  initial  conditions  and  rainfall  inputs. 

Chapter  4  presents  model  predictions  for  a  series  of  observed  rainfall 
events  in  the  Amo  River  basin  in  Tuscany,  Italy. 

Chapter  5  summarizes  the  important  conclusions,  discusses  model 
limitations  and  provides  suggestions  for  future  work. 
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Chapter  2 


MODEL  OF  INFILTRATION 


In  this  Chapter  we  present  the  model  of  infiltration  which  will  be 
used  in  the  rainfall-runoff  model.  The  infiltration  model  is  based  on  the 
kinematic  theory  of  infiltration  which  provides  a  simplified  description  of 
gravity-dominated  moisture  flow.  The  model  incorporates  additional 
approximations  in  the  treatment  of  time-varying  rainfall  rates. 

Kinematic  infiltration  has  been  studied  by  various  researchers 
dealing  with  a  homogeneous,  isotropic  soil  (v.g.,  Smith  and  Hebbert,  1983, 
and  Charbeneau,  1984).  In  Section  2.A  we  review  the  assumptions  of  the 
kinematic  model  of  infiltration  and  the  derivation  of  an  equation  of 
evolution  of  the  depth  of  the  wetting  front  in  a  homogeneous,  inotropic  soil. 
Most  natural  soils,  however,  have  decreasing  permeability  with  depth.  In 
Section  2.B  we  apply  the  kinematic  assumption  to  infiltration  in  a 
vertically  heterogeneous  soil  where  conductivity  decreases  exponentially 
with  depth.  We  consider  both  unsaturated  and  saturated  infiltration  on  a 
hillslope  with  areally  uniform  conditions  of  soil  type,  slope,  and  rainfall 
rate.  Spatial  uniformity  permits  the  use  of  a  one-dimensional  continuity 
equation.  The  analysis  clarifies  the  influence  of  slope,  soil  parameters, 
and  rainfall  rate  on  infiltration  and  the  generation  of  infiltration-excess 
runoff  and  interflow.  In  Section  2.C,  a  conceptual  approximation  to  the 
soil  moisture  profile  is  introduced  that  allows  for  computational  simplicity 
in  the  treatment  of  time-varying  rainfall  rates. 
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2.A  THEORY  OF  THE  KINEMATIC  MODEL  OF  INFILTRATION 


This  section  presents  the  kinematic  theory  of  infiltration  for  one¬ 
dimensional  vertical  flow  in  a  homogeneous,  isotropic  soil.  The  theory  of 
kinematic  waves  was  developed  by  Lighthill  and  Whitham  (1955)  for 
estimating  the  propagation  of  flood  waves.  Wave  motion  follows  both 
kinematic  and  dynamic  laws,  and  its  complete  description  therefore 
includes  both  advective  and  dissipative  terms.  The  kinematic  model  of 
wave  motion  considers  only  the  advective  component  of  flow,  neglecting  the 
second-order  spatial  derivatives  that  describe  dispersive  transport.  In  an 
analogous  fashion,  the  kinematic  model  of  infiltration  relies  on  reducing 
the  moisture  flow  equation  to  its  gravitational  component,  i.e.,  the  pore 
pressure  gradient  component  of  flow  is  neglected. 

The  familiar  Darcy's  equation  for  one-dimensional,  vertical 
unsaturated  flow  (where  z  is  positive  downward),  is: 

qz(0)  =  K(0)  •  (  1  -  )  (2.1) 

where  0  is  the  volumetric  soil  moisture  content; 

qz(0)  is  the  downward  unsaturated  flow; 

K(0)  is  the  unsaturated  hydraulic  conductivity, 

and  vj/(0)  is  the  pore  pressure. 

The  first  term  in  the  parenthesis  represents  the  gradient  of 
positional,  or  gravitational  potential.  The  second  term  in  the  parenthesis 
is  the  gradient  of  pore  pressure.  The  gravitational  flux  component,  K(0), 
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depends  only  on  moisture  content,  and  is  the  advective  term  in  the 
equation.  The  flux  induced  by  the  pressure  gradient  or  capillary  flux, 
-K(0M3\y(0)/3z),  depends  not  only  on  mass  concentration  but  also  on  its 
spatial  gradient,  and  is  the  diffusive  term. 

The  kinematic  model  reduces  Equation  (2.1)  to 

qz(0)  =  K(0)  (2.2) 

The  continuity  equation  for  one-dimensional  vertical  flow  is 

30  3qz(0) 

~T~  +  ^  —  0 

at  3z  (2.3) 


Substitution  of  the  Kinematic  Flow  Equation  (2.2)  in  the  Continuity 
Equation  (2.3)  leads  to 


ae  dK(0)  aa 

dt  +  d0  dz 


(2.4) 


Solving  (2.4)  for  dK(0)/d0  with  0  constant  we  obtain 
dK(0)  99  dz  dz  | 

d0  ”  9t  '  30  ~  dt  1  9  conslam  (2'5 

Provided  that  0(z,t)  is  continuous  throughout  the  z-t  space,  Equation 
(2.5)  states  that  the  mobility  of  a  point  of  constant  0  in  the  soil  moisture 
profile  is  given  by  dK(0)/d0.  Integrating  (2.5)  yields  the  expression  which 
is  the  basis  of  the  method  of  characteristic  solutions  to  (2.3), 
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(2.6) 


z(0)  = 


dK(0) 

d0 


(t-x) 


where  x  is  the  time  of  initiation  of  rainfall 


Most  empirical  models  of  the  K-0  relation  are  either  exponential  or 
power-law  relations  with  the  exponent  greater  than  or  equal  to  3.0.  Thus, 
dK(0)/d0  increases  with  0,  i.e.  the  mobility  is  greater  for  higher  9  values. 
This  leads  to  the  important  conclusion  that  because  moisture  content 
decreases  with  depth  at  the  wetting  front,  a  kinematic  front  steepens  with 
time.  Reversely,  in  a  "drainage  wave"  the  moisture  content  increases  with 
depth,  and  the  profile  is  self-spreading. 

The  steepening  nature  of  the  kinematic  wetting  front  eventually 
leads  to  a  discontinuity  in  the  moisture  profile.  Charbeneau  (1984), 
following  Whitham  (1974),  derives  an  expression  for  the  depth  at  which 
the  moisture  gradient  becomes  infinite. 

The  moisture  profile  of  the  true  front  is  continuous  with  depth.  The 
sharpening  tendency  of  the  kinematic  component  of  flow  is  counteracted  by 
the  spreading  effect  of  the  diffusive  component.  The  equilibrium  shape  of 
the  front  depends  on  the  importance  of  diffusion  in  the  particular  soil. 
Figure  2.1  represents  a  true  wetting  front  and  its  kinematic 
approximation. 

Equation  (2.4)  does  not  apply  to  a  discontinuous  profile  since  neither 
dQ/dt  nor  90/9z  are  continuous  at  the  front.  Instead,  we  will  derive  an 
expression  for  the  advance  of  the  kinematic  front  based  on  continuity 
considerations.  From  continuity,  the  sharp  kinematic  wave  must  have  the 
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FIGURE  2.1:  Wetting  fronts  as  defined  by  the  full 
model  and  the  kinematic  model  of  infiltration. 
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same  moisture  volume  as  the  smooth  wave  which  it  approximates. 
Therefore,  the  depth  of  the  kinematic  front,  Zf,  must  satisfy 

rh. 

0  •  dz  =  (z2- Zf)  •  02+ (Zf- Z!>  •  (2.7) 


Integrating  the  Continuity  Equation  (2.3)  over  depth,  we  obtain 


d  9  r*2 

dtl, 0  dz  +  ad,,  q*(e)  dz  =  0 


(2.8) 


Substituting  the  kinematic  flow  Equation  (2.2), 


d  f2* 

*1  9dz 


+  K(0)  I  =0 


(2.9) 


Replacing  the  integral  in  (2.9)  by  (2.7)  and  differentiating  with 
respect  to  time  we  obtain 


dZf  _  K(0j)  -  K(02) 

dt  0j  -  02 


(2.10) 


Integration  of  (2.10)  yields 


fK(0!)-K(02)A  /  x 

M-'O-er  J-(t-T) 


(2.11) 


where  x  is  the  time  of  initiation  of  rainfall. 
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Note  that  given  the  equal  volume  of  water  in  the  true  front  and  the 
kinematic  front  (Equation  (2.7)),  the  rate  of  advance  of  the  true  wetting 
front  equals  that  of  the  kinematic  front,  given  in  (2.9).  For  our  purposes, 
the  mean  position  of  the  wetting  front,  not  its  shape,  is  of  concern. 


2.B  KINEMATIC  INFILTRATION  IN  A  VERTICALLY 
HETEROGENEOUS  SLOPED  SOIL 


In  this  section,  we  consider  infiltration  in  an  isotropic  hillslope  soil 
where  saturated  hydraulic  conductivity  decreases  exponentially  with 
depth,  under  areally  uniform  soil  parameters  and  slope.  The  analysis  is 
extended  to  anisotropic  soils  in  Appendix  A.  Rainfall  rate,  R,  is 
considered  to  be  constant  in  space  and  time.  The  coordinate  directions  and 
angles  required  by  our  analysis  are  defined  in  Figure-caption  2.2. 


i)  Soil  parameterization  and  moisture  initialization 


Saturated  hydraulic  conductivity  decreases  exponentially  with  depth 
measured  from  the  surface  in  the  vertical  direction  i.e.,  in  the  z  direction. 


41 


FIGURE  2.2:  Representation  of  the  coordinate  directions  on  a 
hillslope  of  constant  slope,  tan(a).  Axis  p  is  parallel  to  the 
hillslope  surface  and  points  in  the  direction  of  maximum  slope. 
Axis  n  is  normal  to  the  hillslope  surface.  Axis  h  is  the  projection 
of  axis  p  on  the  horizontal  plane.  Axis  z  is  in  the  vertical 
direction.  Infiltration  may  be  at  an  angle  a'  with  respect  to  the 
vertical  direction,  as  represented. 
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(2.12) 


Ks(z)  =  Kq-  e-f  z 

where  Ks(z)  is  the  saturated  hydraulic  conductivity  at  depth  z; 
Ko  is  the  surface  saturated  hydraulic  conductivity; 
and  f  is  a  parameter  with  dimensions  L‘l. 


The  Brooks-Corey  (Brooks  and  Corey,  1964)  parameterization  of 
unsaturated  hydraulic  conductivity  will  be  used, 


(2.13) 


where  K(0)  is  the  hydraulic  conductivity  at  moisture  content  0; 

Ks  is  the  saturated  hydraulic  conductivity; 

0S  is  the  saturated  moisture  content; 

0r  is  the  residual  moisture  content,  defined  as  the  value  below 
which  moisture  cannot  be  extracted  by  capillary  forces;  and 
e  is  a  pore  size  distribution  index. 


A  correlation  normally  exists  between  each  of  the  parameters  0r,  0S, 
and  e  and  the  saturated  hydraulic  conductivity.  For  saturated  conductivity 
varying  with  depth,  these  parameters  should  correspondingly  be  functions 
of  depth.  Not  knowing  those  functions,  and  for  simplicity,  we  consider 
parameters  0r,  0S,  and  e  to  be  constant  with  depth. 

Substituting  (2.12)  into  (2.13),  we  obtain 


K(0,z)  =  Kq-  e  fz- 


0  -0,  Y 
0-0r 


(2.14) 


It  is  assumed  that  prior  to  the  initiation  of  rainfall  the  soil  had  a 
constant  recharge  rate,  Rj.  The  initial  moisture  profile  will  be  computed 
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from  this  infiltration  rate.  Estimation  of  Rj  will  be  dealt  with  in  Section 
3.C. 


ii)  Unsaturated  infiltration 


If  the  rainfall  rate  is  less  than  the  surface  saturated  conductivity 
(i.e.,  R<K0),  then  infiltration  occurs  initially  under  unsaturated 

conditions.  For  the  kinematic  approximation,  flow  is  in  the  direction  of  the 
gravitational  potential  gradient  i.e.,  it  is  vertical.  Therefore,  we  have  for 
unsaturated  infiltration,  qz(9,z)  =  Kj.(9,z)  =  R  (for  z  <  Zf).  However,  soil 

parameterization  (2.12)  implies  that  any  rainfall  rate  R  below  the  surface 
saturated  conductivity  has  a  corresponding  critical  depth,  Z*(R),  in  the  soil 
profile  for  which  the  saturated  conductivity  equals  R,  i.e., 

z<  Z*  ;  K,(z)  >  R 

z  =  Z*  ;  K,(z)  =  R 

z  >  Z*  ;  Kg(z)  <  R 


Letting  z=Z*(R)  in  (2.12)  and  solving  for  Z*(R),  we  obtain 

Z*(R)  =  j  •  ln(-^  J  (2.15) 
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The  wetting  front  is  unsaturated  for  Zf  <  Z*(R)  and  saturated  for  Zf  > 

Z*(R).  Saturated  infiltration  is  the  subject  of  section  iii).  Substituting 
K(0,z)  =  R  into  (2.14)  and  solving  for  0(R,z),  we  obtain 

0(R,z)  =  j  £ '  eT z  ‘  (0s“9r)  +  ,  for  z  <  Zf  with  Zf  <  Z*(R)  (2.16) 


Equation  (2.16)  shows  that  the  moisture  content  above  the  wetting 
front  increases  exponentially  with  depth.  Figure  2.3  represents  moisture 
profiles  for  different  rainfall  rates. 

Since  unsaturated  flow  is  vertical,  the  one-dimensional  flow 
continuity  equation  applies.  The  derivation  of  the  expression  of  evolution 
for  Zf  is  analogous  to  that  for  unsaturated  flow  in  a  homogeneous  soil 

presented  in  2.A.  Integrating  (2.3)  in  z  , 

~  0  dz  +  qz  (0,z)  |  =0 

atJZl  Zj 

or 

ijy-dz  +  R-  Rj  =  0  (2.17) 


Through  substitution  of  (2.16),  the  integral  of  moisture  in  (2.17)  is 
given  by 
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FIGURE  2.3:  Moisture  profiles  for  various  rainfall 
rates  lower  than  the  soil  surface  saturated 
conductivity,  for  a  soil  where  hydraulic 
conductivity  decreases  exponentially  with  depth. 
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0(R,z)-dz  + 


0(Ri,z)-dz 


l 


=  j  •(0s-0r)*(j)(e7  ^  ~  el  Z0  +  ©^(Zf-Zj)  + 
E  •(0s-0r){j)-(e  e -  efz0  +  0r-(ZrZl) 


(2.18) 


Differentiation  of  (2.18)  with  respect  to  time  yields  the  first  term  in 

(2.17), 


d  f*2  dZf 

dtj,  dt 


i 

dZf  r  (  R:  \  c 

“((kJ'J  '(0s*0r>ee  f+0r 


dZf 

=  -=L-[0(R,Zf)  -  0(Ri,Zf)l 


(2.19) 


Substituting  (2.19)  into  (2.17)  and  solving  for  dZ/dt  yields 


dZf 

“dt~ 


R-Ri 


0(R,Zf)  -  0(Rj,Zf) 


(2.20) 
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in)  Saturated  infiltration 


The  wetting  front  is  saturated  in  two  cases;  1)  when  the  rainfall  rate 
is  higher  than  the  surface  saturated  conductivity,  or  2)  when  the  wetting 
front  has  penetrated  beyond  the  critical  depth  Z*(R).  In  the  first  case,  the 
entire  wetted  soil  (i.e.,  from  the  surface  to  Zf)  is  saturated.  In  the  second 

case,  as  the  front  reaches  Z*(R),  discharge  can  only  be  less  than  recharge 
from  above,  and  moisture  progressively  accumulates  above  the  wetting 
front.  A  zone  of  perched  saturation  develops  and  grows  upward  from 
Z*(R),  as  well  as  downward  as  the  front  progresses.  Here  we  introduce  a 
new  variable,  Zt,  defined  as  Zf  minus  the  depth  of  the  zone  of  perched 

saturation.  If  there  is  a  zone  of  perched  saturation,  this  definition 
corresponds  to  the  vertical  distance  from  the  top  of  the  zone  of  saturation  to 
the  soil  surface.  If  the  whole  wetted  soil  is  unsaturated  (Zf  <  Z*(R)),  then  Zt 
equals  Zf.  An  equation  of  evolution  for  Zt,  as  well  as  for  Zf,  will  be  derived 
below.  Figure  2.4  represents  the  progression  of  the  wetting  front  in  time, 
showing  the  development  of  a  zone  of  saturation,  and  the  upward  evolution 
of  Zt. 

The  flow  at  the  wetting  front  is  in  the  direction  of  the  gravitational 
potential  gradient,  i.e.  it  is  vertical.  Neglecting  pressure  gradients  we 
have, 


qz(Zf)  =  Ks(Zf)  =  K0-  e"f  Zf  ,  forZf>Z*(R)  •  (2*21) 

Above  the  wetting  front  (z  <  Zf),  the  saturated  conductivity  is  higher 
than  at  the  wetting  front,  but  flow  is  constrained  by  continuity.  Thus,  the 
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FIGURE  2.4:  Evolution  of  the  moisture  profile  over  four  consecutive  time 
steps,  ti  through  t4,  under  constant  rainfall  rate,  R.  At  time  ti  the 
wetting  front  is  unsaturated.  At  time  t2,  Zf  reaches  the  critical  depth 
Z*(R)  for  which  we  have  KS(Z*)  =  R.  From  time  t£  on,  a  zone  of 
saturation  develops  and  grows.  At  time  t4,  the  top  of  the  zone  of 
saturation  reaches  the  soil  surface,  i.  e.  we  have  Zt  =  0. 
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flow  above  the  wetting  front  is  deflected  laterally  and  is  at  an  angle,  a'(z), 
with  the  vertical  direction.  The  flow  deflection  angle  with  respect  to  the 
normal  direction,  i.e.  angle  (a  +  a'(z))  in  Figure  2.5,  is  related  to 
conductivity,  Ks(z),  through  the  law  of  refraction  (see  e.g.  Freeze  and 

Cherry,  1979,  pp  172-3); 

KjCz^  tanfo+a’Cz!)) 

K M  =  tan(a+a'(z2))  ,  for  Zt  £  zh  z2  <  Zf  (2.22) 


The  boundary  condition  at  the  wetting  front  is 
Es(Zf)  =  Ko.e-fZf,  a'(Zf)  =  0 

Substitution  of  (2.23)  into  (2.22)  yields 

Ks(z)  tan(a+a'(z)) 

Ks(Zf)  =  tan(a) 

Substituting  (2.12)  and  solving  (2.24)  for  a  (z), 
a’(z)  =  tan-1[  ef(Zrz),tan(a)  ]  -  a 


(2.23) 


(2.24) 


(2.25) 


Equation  (2.25)  shows  that  the  angle  of  flow  with  respect  to  the 
vertical  direction  increases  upward  from  Zf  to  the  top  of  the  perched 
saturation  zone  (depth  Zt).  Figure  2.6  represents  angle  a’(z)  over  depth 
within  the  zone  of  saturation,  given  f,  Zt  and  three  different  values  of  Zf. 
Since  a'(z)  decreases  with  depth,  the  flow  lines  within  the  zone  of 
saturation  are  curved,  as  represented  schematically  in  Figure  2.7. 
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FIGURE  2.5;  Representation  of  the  flow  vector,  qs, 

and  the  component  of  flow  in  the  normal  direction, 
qn.  The  angle  of  flow  with  respect  to  the  normal 

direction  is  (a+a'(z)). 
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)  over  the  saturated  depth,  with 
and  a  =  45°,  for  three  wetting 
Zf  =  600  mm  and  Zf  =  1000  mm. 


FIGURE  2.7:  Flow  lines  in  the  unsaturated  and  saturated  zones 
of  the  wetted  soil.  Flow  is  vertical  in  the  unsaturated  zone.  In 
the  saturated  zone,  decreasing  permeability  with  depth  causes 
deflection  of  flow.  Flow  deflection  decreases  with  depth  and 
flow  is  vertical  at  the  wetting  front.  For  large  (Zf  Zt)  or  large  f 
values,  flow  near  the  top  of  the  saturation  zone  (i.  e.  at  depth  Zt) 
is  nearly  parallel  to  terrain  slope. 
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Figure  2.8  represents  a'(Zt)  as  a  function  of  terrain  inclination  angle  a,  for 
two  different  values  of  f  and  given  values  of  Zf  and  Zt.  We  see  that  for  large 
f  values,  the  angle  approaches  a'(z)  =  90°  -  a,  that  is,  flow  becomes  closely 
parallel  to  the  terrain  surface. 

At  a  depth  z  within  the  zone  of  saturation,  flow  has  both  a  vertical 
and  a  horizontal  component,  q^z)  and  qy^z),  such  that, 


tan(a'(z))  = 


Qhk) 

qz(z) 


,  for  Zt  £  z  £  Zf 


(2.26) 


Given  that  the  soil  is  isotropic,  the  ratio  between  the  horizontal  and 
vertical  components  of  flow  equals  the  ratio  between  the  hydraulic 
gradients  in  these  two  directions,  Jh(z)  and  Jz(z).  Thus,  from  (2.26), 


.  .  ..  Jh(z) 
tan(a'(z))  =  -7 -7 
Jz(z) 


,  for  Zt  5  z  <  Zf 


(2.27) 


The  vertical  gradient  of  hydraulic  potential  is 


J2(z)h- 


9d>(z,h) 

dz 


d_ 

dz 


e(z,h)  + 


P(z)  ' 
Y  . 


1  dP(z) 
Y  dz 


(2.28) 


where  0(z,h)  is  the  total  gravitational  potential  at  location  (z,h); 

e(z,h)  is  the  elevation  potential  at  location  (z,h); 

and  P(z)/y  is  the  pressure  potential  induced  by  gravity  at  depth  z. 


The  horizontal  gradient  of  hydraulic  potential  is 


4(z)*“ 


3<l>(z,h) 

ah 


a 

ah 


e(z,h)  + 


P(z)  ' 
Y  . 


1  aP(z) 

Y  ah 
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FIGURE Angie  of  flow  at  the  top  of  the  zone  of 
saturation,  ct'(Zt),  as  a  function  of  terrain  inclination 

angle  a,  with  Zt  =  o  mm,  Zf  =  1000  mm,  and  (a)  f  =  io*2  mm'^; 
and  (b)  f  =  10'2  mm'1. 
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1  9z  dP(z) 
Y  9h  dz 


t  M1  dP(z) 
tan(q)-- — ^ — 
Y  dz 


(2.29) 


Thus,  we  have 
Jh(z)  =  (1-Jz(z))  •  tan(a) 

Solving  (2.27)  for  Jh(z), 
Jh(z)  =  tan(q'(  J)  •  Jz(z) 


(2.30) 


(2.31) 


Equating  (2.30)  and  (2.31)  we  obtain  J2(z)  and  the  vertical  component 
of  flow; 


Jz(z)  = 


tan(a) 

tan(a)  +  tan(a'(z))  ; 


qz(z)  =  Jz(z)-K0e_fz 


(2.32) 


Substituting  (2.32)  into  (2.31)  we  obtain  Jh(z)  and  the  horizontal 
component  of  flow; 


Jh(z)  = 


tan(q)-tan(a'(z)) 
tan(q)  +  tan(q'(z))  ’ 


qji(z)  =  Jh(z)-K0-e“f'2 


(2.33) 


Figures  2.9  and  2.10  represent  the  vertical  and  horizontal  gradients 
over  depth  within  the  zone  of  saturation,  for  given  f,  Zt  and  three  different 
values  of  Zf.  Figures  2.11  and  2.12  represent  Jz(Zt)  and  Jh(Zt)  as  a  function 
of  terrain  inclination  q,  for  given  f,  Zt  and  Zf.  We  see  that  Jh(Zt)  has  a 
maximum  at  q  =  45°.  This  indicates  that  maximum  lateral  flow  occurs 
with  q  =  45°. 
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FIGURE  2.9;  Profile  of  Jz(z)  over  the  saturated  depth,  with  Zt 

=  o  mm,  f  =  10*3  mm*1,  and  a  =  45°,  for  three  wetting  front 
depths,  Zf  =  200  mm,  Zf  =  600  mm  and  Zf  =  1000  mm. 


Jh 


EKrUBE  2.X0«  Profile  of  J^(z)  over  the  saturated  depth,  with 

Zt  =  o  mm,  f  =  10'3  mm-1,  and  a  =  45°,  for  three  wetting 
front  depths,  Zf  =  200  mm,  Zf  =  600  mm  and  Zf  =  1000  mm. 
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■FIGURE  2.1.1 :  V erticai  gradient  of  hydraulic  potential  at  the 
top  of  the  zone  of  saturation,  Jz(Zt),  as  a  function  of 

terrain  inclination  angle  a,  with  Zt  =  0  mm,  Zf  =  1000  mm, 
and  f  =  io*3  mm*1. 
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FIGURE  2.12:  Horizontal  gradient  of  hydraulic  potential  at 
the  top  of  the  zone  of  saturation,  Jh(Zt),  as  a  function  of 

terrain  inclination  angle  a,  with  Zt  =  0  mm,  Zf  =  1000  mm, 
and  f  =  10'3  mm'V 
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The  flow  continuity  equation  must  now  be  written  in  two 
dimensions.  As  we  will  see,  it  is  more  convenient  to  write  the  continuity 
equation  for  the  orthogonal  components  of  flow  in  the  n  and  p  directions 
than  for  the  z  and  h  directions.  The  continuity  equation  is 


39  3qn(z)  dqp(z)  ^ 

dt  +  dn  +  dp 


(2.34) 


Integration  of  (2.34)  over  2  (with  Zt  <  z^  <  Zf  and  Z2  >  Zf)  leads  to 
d  f22  d  f22  d  fz2 

dfj2,  9dz  +  tej2l  q"(z>dz  +  aFJ2,  q»(zHz  =  °  (2-35> 


Applying  the  chain  rule, 


d  f^2  dz  d  h  ..  .  dz  df^  A 

dtl  e  dz  + q"(z)dz+3?-alJ2  q»(z)dz=0 


(2.36) 


As  we  have  dz/dp  =  0,  the  last  term  in  (2.36)  is  zero,  and  the  equation 
reduces  to  the  one-dimensional  form  (with  dz/dn  =  l/cos(a)), 

iJ2/dz+^k)'q"(z)C=o  (2-37> 

Within  the  zone  of  saturation,  we  have  0(z)  =  0S  and  d0/dz  =  0.  Since 
2  is  constant  along  direction  p,  we  have  dqp(z)/dp  =  0.  Therefore,  the 
Continuity  Equation  (2.34)  gives  dqn(z)/dn  =  0  within  the  saturated  zone. 
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Thus,  qn(z)  is  constant  in  the  saturated  zone,  i.e.,  qn(z)  =  qn(Zf)  for 
Zt  <  z  S  Zf.  Substituting  qn(Zf)  for  q^zj)  in  (2.37),  we  obtain 


8*dz 


+ 


1 

cos(a) 


•[qn(z2)  -  qn(Zf)]  =  0 


(2.38) 


From  geometrical  considerations,  the  normal  component  of  vertical 
flow  is  qn(z)  =  cos(a).qz(z),  and  (2.38)  can  be  written  as 


■^J  0-dz  +  [c^  (z2)  -  q2(Zf)]  =  0 
—J"2  0dz  +  [Ri-K0e'fZf]  =0 


(2.39) 


where  Rj  is  the  initial  infiltration  rate.  The  procedure  for  estimating  Rj  is 
described  in  Section  3.D. 

Through  substitution  of  Equation  (2.16)  for  8(Rj,z),  the  integral  of 
moisture  in  (2.39)  is 


Differentiating  (2.40)  with  respect  to  time, 
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1 


eH  d2^f 

7  t  I 

{  ^-TT 

e3- 

iKo; 

•(0s-0r)eT  Zf-  0r 

dZf 

=  =-(0.-0(Ri*)) 


(2.41) 


Substituting  (2.41)  into  (2.39)  and  solving  for  dZf/dt,  we  obtain 
dZf  K0,e-f'Zf  -  Rj 

dt  0S  -  0(Rj,Zf)  (2.42) 

Derivation  of  the  equation  of  evolution  for  Zt  is  done  in  entirely 
analogous  manner  to  the  derivation  of  the  equation  of  evolution  for  Zf. 

Integration  of  Equation  (2.34)  over  z  again  leads  to  Equation  (2.37),  but  now 
with  zi  and  z2  defined  by  <  Zt  <  z2  <  Zf.  We  saw  that  qn(z)  is  constant  for 
Zf  £  z  <  Zf.  Substituting  qn(Zf)  for  qn(z2)  in  (2.37),  we  obtain 

9d2  +  ^Ttq"(Zf) '  q"(z‘)1 = 0  <2-43> 


Flow  is  vertical  at  z1  (since  it  is  unsaturated  and  we  are  neglecting 
pore  pressure);  and  is  also  vertical  at  the  wetting  front.  Therefore,  the 
component  of  flow  in  the  n  direction  is  obtained  geometrically  from  the 
flow  as  qn(z)  =  cos(a).qz(z),  and  (2.43)  becomes 

jyj2  9‘dz  +  [q*(Zf)  -  qz(z1)]  =  0 
or 
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(2.44) 


0-dz  +[Ko-e"fZf-R]  =  0 


The  integral  of  moisture  in  (2.44)  is 


dz  =  0s-(z2  -  Zt)  + 


f 

(0S— 0r)'e  e  z  +  0r 


•  dz 


=  0S  (Z2  -  Zt)  + 


f 

-  e” z')  +  0r-(Zt-z1)(2.45) 


Differentiating  (2.45)  with  respect  to  time, 


l 


r  f  i 

i, 9dz = tt 

.(koJ 

(0s-0r)eT'z  +  0r  -  0g 

dz, 

=  -^-(0(R,Zt)  -  0a) 


(2.46) 


Substituting  (2.46)  into  (2.44)  and  solving  for  dZt/dt  yields 

dZt  Ko-e~fZf-R 
“dT=  0S  —  0(R,Zt) 


Since  we  have  R  >  KQ.e'^-2^  the  time  differential  of  Zt  is  negative,  i.e., 
Zt  approaches  the  terrain  surface.  Eventually,  Zt  will  reach  the  surface  (Zt 
=  0),  and  from  then  on  we  must  have 


dZ. 

Zt  =  0 ;  -dT  =  0 


(2.48) 
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For  Zt  =  0,  compatibility  of  (2.47)  with  (2.48)  requires  that 
R  =  Ko.e'f  zf,  with  R  now  representing  the  rate  of  infiltration.  Therefore, 
when  Zt  reaches  zero,  the  rate  of  infiltration  drops  from  the  rainfall  rate  to 

the  saturated  conductivity  at  the  wetting  front.  Designating  the  rate  of 
infiltration  by  Rinf  and  the  rate  of  runoff  by  Rn 


^inf  - 


(Before  ponding  (Zt  >  0):  R 

After  ponding  (Zt  =  0):  K<).e'f'zf 


Rr  -  R  -  Rinf 


(2.49) 


2.C  VARIABLE  RAINFALL  RATES 


During  a  rainstorm,  the  rainfall  rate  at  any  location  varies  in  time. 
This  variation  can  be  approximated  by  a  succession  of  short  time  intervals 
of  steady  rainfall  rate.  At  the  start  of  each  time  interval,  a  new  wave  of 
moisture  leaves  the  surface  and  propagates  downward  in  the  soil  profile. 
The  vertical  moisture  profile  can  therefore  be  conceptualized  as  a 
succession  of  waves.  Figures  2.13  and  2.14  represent  complex  moisture 
profiles  for  a  soil  where  conductivity  decreases  exponentially  with  depth. 
If  a  rainfall  rate  Rj  is  succeeded  by  a  higher  rainfall  rate  R2,  the  moisture 

profile  will  initially  resemble  the  schematic  representation  of  Figure  2.13. 
If  instead  we  have  R2  being  less  than  Ri,  then  the  moisture  profile  will 
resemble  that  in  Figure  2.14. 
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FIGURE  2.13:  Moisture  profile  after  a  step  increase 
in  rainfall  rate  for  a  soil  where  conductivity 
decreases  exponentially  with  depth. 


FIGURE  2.14:  Moisture  profile  after  a  step  decrease 
in  rainfall  rate  for  a  soil  where  conductivity 
decreases  exponentially  with  depth. 
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Modeling  of  this  process  of  successive  moisture  waves,  propagating 
at  different  velocities,  overruning  each  other  or  spreading  apart,  would 
require  keeping  track  of  multiple  state  variables,  each  representing  the 
depth  Zfj  of  each  wave  i.  Charbeneau  (1984)  presents  solutions  for  the  soil 

moisture  profile  for  a  succession  of  time  intervals  of  different  rainfall 
rates.  However,  both  memory  and  computational  requirements  to  carry 
out  such  computations  would  be  prohibitive.  For  simplicity,  a  conceptual 
approximation  is  introduced  to  the  vertical  moisture  profile.  The  complex 
profile  of  succeeding  waves  is  approximated  by  a  single-wave  profile  that 
contains  the  same  volume  of  moisture. 

The  single-wave  approximation  requires  introducing  a  third  state 
variable,  representing  the  uniform  rate  of  vertical  flow  that  corresponds  to 
the  single  wave  approximation.  The  new  variable  is  designated  Re,  and 
can  be  conceptualized  as  representing  an  “ equivalent "  steady  rainfall  rate 
in  the  sense  that  the  kinematic  wave  that  corresponds  to  rate  Re  contains 
the  same  volume  as  the  true,  complex  profile.  (Note  that  Re  generally  does 

not  coincide  with  the  average  rainfall  rate,  given  the  non-linearity  of  the 
K-0  relation.) 

We  can  derive  an  expression  for  Re  from  the  moisture  storage  in  the 
soil.  Re  is  a  function  of  the  volume  of  moisture  stored  in  the  soil  above 
depth  Zt  i.  e.,  in  the  unsaturated  zone  of  the  wetted  soil.  The  moisture 
volume  per  unit  area  (dimension  L)  above  the  wetting  front  is  Mt  (for 
"total"  moisture).  The  moisture  volume  per  unit  area  above  Zt  is  Mu  (for 
"unsaturated"  moisture).  Mu  is  obtained  from  Mt  by  subtraction  of  the 
moisture  volume  in  the  saturated  zone  between  depths  Zf  and  Zt.  (Note 
that  if  the  front  is  unsaturated,  we  have  Zt  =  Zf  and  Mu  =  Mt.)  Thus, 
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(2.50) 


(2.51) 

(2.52) 


By  definition  of  Re,  the  single-wave  approximation  0(Re,z)  contains 
the  same  volume  of  moisture  as  the  complex  true  moisture  profile  0(z). 
Therefore,  we  may  replace  0(Re,z)  for  0(z)  in  the  above  integrals. 

Expression  (2.51)  becomes: 


•  dz  =  Mu 


(2.53) 


This  integral  can  be  solved  for  Re  analytically  upon  substitution  of 
the  expression  for  0(Re,z).  The  expression  for  0(R,z)  is  given  in  (2.16). 
Letting  R  =  Re,  and  substituting  (2.16)  into  (2.53)  we  obtain 

-T  i 

■dz  =  Mu 

(2.54) 


Jo 


1  Y  R  f 

(k^J  -(0s-0r)-^'Z  + 


Solving  for  Re, 


Re  =  Ko- 


Mu  -  0r-Zt 


(0s-0r){y)(efZ'-I) 


(2.55) 
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Note  that  if  Zt  is  at  the  surface  i.  e.,  if  the  full  depth  of  wetted  soil  is 
saturated,  there  is  no  unsaturated  moisture  profile  to  approximate  (Mu  = 
0).  In  this  case,  we  define  Re  =  R.  The  expressions  derived  above  for  Re 
apply  for  Zt  >0. 
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Chapter  3 


THE  RAINFALL-RUNOFF  MODEL  AT  THE  BASIN  SCALE 


In  this  Chapter  we  present  a  rainfall-runoff  model  to  be  used  for 
flood  prediction  at  the  basin  scale.  The  rainfall-runoff  model  considers  the 
soils  in  the  basin  to  be  represented  by  the  parameterization  of  hydraulic 
conductivity  given  in  Equation  (2.12),  and  uses  the  kinematic  model  of 
infiltration  presented  in  Chapter  2. 

In  Chapter  2  we  considered  infiltration  under  idealized  spatial 
uniformity  of  all  parameters  and  variables.  Over  a  catchment  as  well  as  at 
the  hillslope  scale,  however,  soil  properties,  slope,  moisture  conditions, 
and  rainfall  inputs  vary  areally.  In  Chapter  2  we  saw  that  the  kinematic 
model  of  infiltration  predicts  lateral  subsurface  saturated  flows  in  a  soil 
where  conductivity  decreases  with  depth.  Under  the  condition  of  spatial 
uniformity,  kinematic  infiltration  was  described  by  the  one-dimensional 
form  of  the  continuity  equation.  Under  spatially  heterogeneous  conditions, 
however,  infiltration  must  be  considered  in  three  dimensions.  Given  that 
full  solution  of  three-dimensional  differential  equations  would  be 
computationally  prohibitive  at  the  basin  scale,  the  model  uses  spatial 
discretization  into  areally  uniform  elements  to  reduce  three-dimensional 
differential  equations  to  their  one-dimensional  form.  Thus,  the  infiltration 
equations  developed  in  Chapter  2  apply  within  each  element.  Spatial 
discretization  is  described  in  Section  3. A.  Lateral  flows  between  elements 
is  accounted  for  through  a  computational  scheme  that  allows  simplified 
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element  coupling.  The  element  coupling  procedure  is  described  in  Section 
3.B. 

The  computational  procedure  also  involves  temporal  discretization. 
At  each  time  step,  the  state  of  moisture  of  each  element  is  updated  and  the 
rate  of  infiltration-excess,  saturation-excess  and  return  flow  runoff  is 
computed. 

To  determine  saturation-excess  runoff,  the  model  considers  an 
estimated  initial  water-table.  The  procedure  for  water-table  estimation 
and  for  moisture  initialization  above  the  water-table  is  described  in  Section 
3.D.  Rainfall  over  elements  where  the  water-table  is  at  the  soil  surface 
(saturated  elements)  becomes  saturation-excess  runoff.  Infiltration- 
excess,  saturation-excess  and  return  flow  runoff  contribute  to  streamflow, 
and  are  routed  over  the  hillslope  and  along  stream-channels  to  the  basin 
outlet.  Section  3.E  describes  the  routing  procedure. 


3A  SPATIAL  DISCRETIZATION  AND  COMPUTATIONS  IN  ONE 
ELEMENT 


Basin  elements  are  taken  as  the  pixels  of  the  D.E.M.  grid.  This  is 
done  for  various  reasons  which  are  summarized  next.  First,  the  use  of 
elements  in  a  regular  grid  offers  considerable  computational  advantage. 
Element  parameters  and  variable  values  can  be  stored  in  a  rectangular 
matrix,  and  each  element  is  easily  referenced  by  its  row  and  column 
numbers,  facilitating  both  compact  data  storage  and  rapid  data  extraction; 
element  boundaries  are  rectangular  and  neighboring  elements  are  easily 
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identified;  computation  algorithms  are  greatly  simplified.  Second,  the 
model  emphasizes  terrain  morphology,  and  the  variables  describing 
topography  can  be  extracted  most  simply  and  without  bias  from  a  regular 
matrix  of  terrain  elevations,  as  opposed  to  through  the  extraction  of 
contour  lines  (Tarboton,  1989).  Third,  algorithms  for  extraction  of  river- 
channel  networks  from  regular  D.E.M.s  have  been  developed  and  proven 
to  compare  well  with  mapped  data,  and  the  use  of  the  D.E.M.  grid  for  basin 
discretization  allows  straightforward  incorporation  of  extracted  channel 
networks.  Finally,  D.E.M.  rectangular  grids  are  available  for  great  part  of 
the  U.S.  territory  (at  30m  x  30m  and  60m  x  90m  grids)  and  many  other 
parts  of  the  World.  The  available  resolution  is  in  general  superior  to  that 
of  soil  maps  (from  which  soil  parameters  are  estimated)  and  to  that  of 
radar  measurements  of  rainfall  rate. 

Each  element  of  the  model  consists  of  a  soil  column  with 
rectangular  horizontal  cross-sectional  area  and  is  characterized  by  its  soil 
properties,  maximum  slope  (tan(a))  and  direction  of  maximum  slope.  Soil 
properties  may  vary  between  elements,  and  in  general  can  be  obtained 
from  digitized  soil  maps. 

It  is  assumed  that  soil  properties,  slope,  rainfall  rate,  and  the  initial 
moisture  profile  are  areally  uniform  within  an  element.  The  equations 
derived  in  the  previous  section  for  a  hillslope  of  spatially  uniform 
characteristics  therefore  apply  within  an  element.  Spatial  variability  of 
soil  type,  terrain  slope  and  rainfall  rate  exists  among  elements,  and 
lateral  subsurface  discharge  of  moisture  from  one  element  to  the  next 
differs  among  elements.  Inputs  and  outputs  of  moisture  into  and  from  an 
element  are  not  equal,  and  lateral  flows  affect  the  moisture  storage  in  any 
element.  Therefore,  even  having  reduced  infiltration  in  each  element  to  its 
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one-dimensional  form,  we  must  account  for  the  balance  of  lateral  moisture 
inputs  and  outputs  when  updating  the  moisture  storage  in  the  element. 
In  Section  3.B,  we  present  a  computational  procedure  that  allows  for 
element  coupling  through  lateral  flows. 

Complete  description  of  the  state  ox  moisture  in  one  element  requires 
three  state  variables.  Mt,  Zf  and  Zt,  defined  in  Chapter  2.  All  other 

variables  required  to  write  the  equations  of  evolution  (i.  e.  the  differential 
equations  in  time)  for  the  state  variables  are  obtained  from  them.  The 
equations  of  evolution  for  Zf  and  Zt  were  derived  in  Section  2.B.  The 
equation  of  evolution  for  the  moisture  storage,  Mt,  is  obtained  from 

straightforward  mass  balance  considerations.  The  input  flows  are 
infiltrated  rainfall,  and  the  lateral  flows  from  contributing  elements,  Q;n. 
The  output  is  the  lateral  flow,  Qout,  from  the  element  of  interest  to  a 
neighboring  element. 

As  we  saw  in  2.B,  lateral  flows  are  generated  only  from  the 
saturated  zone.  Thus,  the  general  expression  for  the  volumetric  discharge 
from  an  element  to  the  element  situated  downslope  from  it  is 

j r2* 

Qout  =  w-  qh(z)'dz  (31) 

where  W  represents  the  width  of  flow,  i.e.,  the  width  perpendicular 
to  the  flow  direction.  For  flow  in  the  x  or  y  directions,  W  is  Ay  or  Ax  , 
respectively;  for  flow  at  an  angle  with  the  x  and  y  directions,  W  is 

(Ax-Ay )/y  Ay2+Ay2. 

Through  substitution  of  Equation  (2.33)  for  q^Cz)  in  the  above 
integral, 
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„  f2^  tan(a)-tan(a’(z))  Tr  _f.z 

^  W  tan(a)  +  tan(a’(z))  K°  dz 


The  above  integral  equals  the  following  (see  Appendix  B), 
Qout=  W-Ko-cos(a)-sin(a)-e‘f'Zf|  Y'(e'f'z‘-  e'fZ0  -  (Zf-Z^ 


(3.3) 


Below  we  summarize  the  equations  of  evolution  for  the  three  state 
variables,  and  the  expressions  through  which  the  other  variables  involved 
in  those  equations  are  obtained  from  the  state  variables.  In  the  equations 
below,  A  designates  element  area.  Estimation  of  Qjn  is  achieved  through 

the  element  coupling  scheme  which  is  explained  in  Section  3.B. 


dZ i 
dt 


dZt 


dZf  Qin  Q( 

~ 0(Ri,Zf)  +  Min(R,K0)  +  --- 

< 

unsaturated: 

Re-Ri 

0(Re,Zf)  -  0(Ri,Zf) 

< 

saturated: 

Ko'e~f  Zf  -  Rj 

0S  -  0(Rj,Zf) 

r 

unsaturated: 

dZf 

dt 

< 

Ko-e“f  Zf-  Re 

saturated: 


0S  -  0(R„Zf) 


where 
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Re  =  Ko- 


r  mu- 

8,-A  V 

(0,-0r)  | 

ri] 

IfJ 

,  i.7 
(et  Zt-  1) 

Mu  =  Mt-es-(ZrZt) 

0(R,z)  =  E  •  e7'z  •  (0s-er)  +  0r 

Qout=  W-Ko-cos(a)-sin(a)-e-fZf  [  j-(e'f'Zt  -  e’fz0  -  (Zf-Z^ 


33  ELEMENT  COUPLING 


When  subsurface  flow  has  a  lateral  component,  the  state  of  elements 
located  upstream  and  downstream  mutually  affect  each  other.  Thus,  the 
equations  of  evolution  of  the  three  state  variables  are  coupled  to  all  other 
equations  of  the  elements  that  contribute  subsurface  flows  to  the  same 
stream  element.  Simultaneous  integration  of  the  large  number  of  coupled 
differential  equations  would  be  exceedingly  numerically  intensive  given 
the  objectives  of  this  work,  but  if  the  lateral  component  of  flow  is 
significant,  computations  should  incorporate  element  interactions,  even  if 
in  an  approximate  way.  We  present  a  computational  procedure  that 
allows  for  element  coupling  through  lateral  flows  in  a  simplified  manner. 
The  computational  technique  treats  each  element  at  a  time  but  processes 
elements  in  an  order  that  still  allows  for  coupling  of  elements  through 
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lateral  flows.  Thus,  just  three  differential  equations,  corresponding  to  the 
state  variables  of  one  element,  need  be  integrated  jointly. 

Computations  start  with  the  uphill  elements  and  work  their  way 
down  to  the  elements  at  the  bottom  of  the  hillslope.  This  order  of 
computations  guarantees  that  when  the  state  variables  are  updated  for  any 
element,  the  outflow  from  each  element  that  discharges  into  it  has  already 
been  computed.  The  computational  procedure  relies  on  time  as  well  as 
spatial  discretization.  Within  a  time  interval,  the  outflow  from  one 
element  is  computed  by  integration  of  the  lateral  flows  in  time.  However, 
the  average  outflow  rate  (total  discharge  divided  by  the  time  i nterval 
considered)  is  used  in  updating  the  state  variables  in  the  element  located 
downstream  from  it,  and  thus  only  one  value  need  be  stored  for  each 
element  in  element  coupling. 

In  updating  the  state  of  moisture  of  each  element,  the  balance 
between  lateral  flows  is  added  to  the  expression  for  the  integral  of  moisture 
above  the  wetting  front,  according  to  the  equations  of  evolution  for  Mt 
provided  in  Section  3. A.  Whenever  the  value  of  Mt  exceeds  the  volume 
corresponding  to  complete  saturation  above  the  wetting  front  (i.e. ,  Zf-0s) 
and  lateral  inputs  exceed  lateral  ouput  from  the  element,  return  flow  is 
generated.  In  other  words,  return  flow  is  produced  when  the  balance  of 
lateral  subsurface  flows  is  positive  and  the  element  is  already  saturated. 
Return  flow  is  added  to  the  surface  runoff  volume  that  is  routed  through 
the  basin. 
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3.C  RAINFALL  INPUTS 


The  rainfall  rate  in  each  element  at  each  time  step  may  be  obtained 
either  from  radar  or  raingage  measurements.  Radar  offers  much  larger 
spatial  and  temporal  resolutions  than  can  be  provided  by  a  network  of 
telemetering  raingages.  The  issue  of  estimation  of  rainfall  rates  from 
radar  reflectivity  measurements  and  its  potential  sources  of  error  are 
reviewed  in  Austin  (1987).  Wyss  et  al.  (1990)  used  radar  estimates  of 
rainfall  as  input  to  a  distributed  rainfall-runoff  model.  When  radar 
measurements  are  not  available,  data  from  a  network  of  telemetering 
raingages  may  be  used.  Various  methods  can  be  employed  to  estimate 
rainfall  rate  distribution  from  the  point  rainfall  measurements.  Perhaps 
the  simplest  of  such  methods  is  to  use  Thiessen  polygons  to  approximate 
rainfall  distribution.  The  method  corresponds  to  approximating  the 
rainfall  rate  in  each  element  by  the  reading  from  the  closest  raingage. 
This  method  will  be  used  in  the  model  application  presented  in  Chapter  4. 


3.D  GROUNDWATER  INITIALIZATION 


To  account  for  the  saturation-excess  mechanism  of  runoff 
generation,  the  event-based  model  must  be  initialized  with  an  estimated 
profile  of  the  water-table.  In  this  section,  we  describe  how  a  computational 
procedure  parallel  to  that  presented  in  Section  2.B  for  the  update  of  the 
near-surface  moisture  profile  during  a  rainstorm,  can  be  used  to  obtain 


77 


the  water-table  position,  and  the  areas  of  soil  saturation,  at  the  time  of 
onset  of  the  rainfall  event. 

The  model  relies  on  a  series  of  assumptions:  (1)  the  rate  of  recharge 
of  the  water  table  is  a  slow-varying  value  with  a  seasonal  cycle;  (2)  the 
water  table  adjusts  to  the  seasonal  fluctuations  of  recharge  rate  rapidly 
enough  that,  at  any  time  between  storms,  it  is  in  a  steady  state  of 
equilibrium  with  the  recharge  rate;  (3)  the  rate  of  water-table  recharge 
between  storms  is  spatially  uniform  over  the  catchment;  (4)  the  aquifer  is 
unconfined  and  covers  the  entire  basin;  and  (5)  groundwater  flow  in  the 
general  zone  of  saturation  can  be  described  as  a  Darcian  flow  system  i.  e., 
macropore  and  pipe  flow  are  negligible.  These  assumptions  are  implicit 
also  in  two  models  that  have  a  similar  intent  to  the  model  presented; 
TOPMODEL  of  Beven  et  al.  (Beven  and  Kirkby,  1979;  Beven  and  Wood,  1983; 
Beven  et  al.,  1985),  and  the  models  of  O'Loughlin  and  Moore  (O'Loughlin, 
1981,  1986;  Moore  et  al.,  1986). 

Freeze  and  Cherry  (1979)  state  that  assumptions  (1)  and  (2)  are 
supported  by  field  observations.  It  seems  reasonable  to  expect  that  these 
assumptions  will  hold  least  strongly  for  low-transmissivity  soils  in  regions 
where  the  range  of  values  of  water-table  recharge  varies  greatly  from  the 
wet  to  the  dry  season,  because  the  drainage  of  such  soils  may  be  too  slow  to 
allow  an  equilibrium  to  be  reached  with  the  low  recharge  rate. 
Assumptions  (1)  and  (2)  should  hold  best  for  the  wet  season,  and  for  the 
purpose  of  flood  forecasting  the  estimation  of  the  water-table  position  in  the 
wet,  not  dry,  season  is  of  interest.  Assumptions  (3)  seems  reasonable 
lacking  information  on  the  spatial  distribution  of  the  water-table  recharge 
rate.  However,  spatial  variations  are  to  be  expected  as  a  of  non-uniformity 
of  rainfall  rates,  vegetation,  etc.  There  appears  to  be  the  possibility  of 
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investigating  the  correlation  between  terrain  elevation,  exposure, 
vegetation  coverage,  etc.,  and  water-table  recharge  rates.  Assumption  (4) 
may  or  may  not  be  met,  but  may  represent  a  reasonable  approximation, 
depending  on  the  true  situation  on  a  particular  basin.  The  computational 
scheme  presented  here  allows  consideration  of  spatially  variable  recharge 
rates,  and  of  unconfined  aquifers  covering  only  part  of  the  basin.  However, 
it  does  not  provide  for  confined  aquifers,  as  it  assumes  a  correspondence 
between  level  of  the  water  table  and  hydraulic  potential.  Assumption  (5)  is 
probably  the  weakest,  as  macropores  and  pipes  are  believed  to  be  present  in 
most  catchments.  This  assumption  is  acceptable  in  basins  where 
macropore  connectedness  and  piping  is  not  well  developed. 

The  equilibrium  assumption  (2)  suggests  that  the  rate  of  recharge 
can  be  estimated  from  the  rate  of  discharge  from  the  basin.  The  rate  of 
recharge  is  estimated  by  dividing  the  volumetric  discharge  from  the  basin 
(measured  at  its  outlet)  by  the  basin  area: 


Ri  = 


Q_ 

Ay 


(3.4) 


where  Rj  is  the  water-table  recharge  rate; 

Q  is  the  observed  baseflow  from  the  basin;  and 
Ay  is  the  basin  area. 


Groundwater  flow  in  the  saturated  zone  below  the  water-table  is 
modelled  in  a  manner  similar  to  the  element  coupling  scheme  presented 
in  3.B.  To  compute  water-table  elevation,  the  saturated  zone  is  initialized 
at  full  saturation  of  the  soil  profile  (the  initial  water-table  coincides  with 
the  terrain  surface),  and  allowed  to  drain  over  time,  under  a  constant  and 
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spatially  uniform  rate  of  recharge,  Rit  computed  from  (3.4).  Unless  soil 
transmissivity  of  the  entire  saturated  soil  column  is  very  low,  the  initial 
discharge  from  the  fully  saturated  soil  is  higher  than  recharge,  and 
groundwater  storage  decreases.  Therefore,  the  water  table  drops  below  the 
soil  surface,  at  least  in  part  of  the  basin.  In  time,  the  decrease  in 
groundwater  storage  leads  to  progressively  lower  discharges.  When 
discharge  rate  reaches  the  value  of  the  recharge  rate,  the  groundwater 
system  remains  stationary,  in  a  state  of  equilibrium.  The  values  obtained 
for  the  depth  of  the  water  table  in  each  element  when  the  steady  state  is 
reached  are  used  in  the  initialization  of  the  event-based  rainfall-runoff 
model. 

Only  one  state  variable ,  the  depth  to  the  water-table,  Zwt ,  is  required 

to  characterize  the  state  of  the  general  zone  of  saturation  in  one  element. 
Like  Zf  and  Zt,  Zwt  is  measured  in  the  vertical  (z)  direction,  positive 

downwards,  from  the  terrain  surface,  and  has  dimension  L. 

The  equation  of  evolution  for  Zwt,  obtained  by  derivation  analogous  to 
that  for  Zt  in  Section  2.B,  is 

Qin  —  Qout  „ 

dZ^  a  R> 

dt  0S  -  e(R,,Zw,)  (3-5) 

where  Qjn  is  the  input  subsurface  flow  into  the  element; 

Qout  is  the  subsurface  discharge  from  the  element; 

and  A  is  the  element  area. 

The  general  expression  for  Qout  is  obtained  through  integration  of 
the  saturated  horizontal  flows  over  the  saturated  depth, 
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Qout  =  wf  qh(z)-dz 

JZwt 


(3.6) 


where  Zjmp  is  the  depth  of  an  impermeable  layer  (the  aquifer 
bottom),  and  W  is  the  width  of  flow  perpendicular  to  the  direction  of 
flow  (the  same  as  in  (3.1)). 


Representing  the  hydraulic  gradient  by  VO,  (3.6)  may  be  written  as 


z*. 


VO-Kne"fz-dz 


=  WVO— ~  (e_f Zwt  -  e  fZimp) 

=  W-VO— e"^  ,  for  Z^p  large 


(3.7) 


Note  that,  since  cpnductivity  decreases  exponentially  with  depth, 
knowledge  of  the  aquifer  depth,  Z\mp,  is  not  required  for  computation  of 
Qout’  provided  that  the  soil  is  relatively  deep  (Zjmp  large). 

Given  that  we  are  under  the  assumption  of  the  aquifer  being 
unconfined,  the  hydraulic  gradient,  VO,  may  be  represented  by  the  slope  of 
the  water-table  surface.  We  approximate  water-table  slope  in  element  i  by 
the  forward  difference  in  water  table  elevations  in  the  direction  of  highest 
water-table  slope, 

V  O  =  tan(a)  -  -Zw‘ ^  ~  ‘  (3.8) 
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where  Zwt  I  j  and  Zwt  I  j  represent  the  depth  of  the  water  table  for  two 
elements  i  and  j,  with  element  j  situated  downstream  from  element 
i. 


Similarly  to  the  infiltration  model,  element  coupling  requires  that 
the  order  of  computations  start  with  uphill  elements  and  work  its  way 
down  to  the  elements  at  the  bottom  of  the  hillslope.  When  groundwater 
storage  is  updated  for  each  element,  Qjn  is  obtained  through  the  addition  of 
the  values  of  discharge,  Qout,  from  all  the  elements  that  contribute  flows 
directly  to  the  element  of  interest,  which  given  the  order  of  computations 
have  already  been  computed. 

An  issue  that  must  be  noted  is  that  the  results  depend  on  the 
boundary  condition  given  by  the  hydraulic  head  (elevation)  of  the  water 
surface  in  the  stream  located  at  the  hillslope  bottom.  The  deeper  the 
stream  channel  and  the  lower  the  stage,  the  deeper  we  may  expect  the 
predicted  water  table  to  be  in  the  vicinity  of  the  channel,  since  the 
hydraulic  gradient  will  be  higher  in  this  location.  Data  on  streamwater- 
surface  elevation  along  the  stream  network  is  generally  not  available.  In 
fact,  in  the  general  situation  no  information  is  available  on  channel  depth, 
and  only  a  few  sparse  streamgages  may  be  expected  to  provide  stage 
readings.  In  this  situation,  we  approximate  streamwater-surface 
elevation  by  terrain  elevation  i.  e.,  we  assume  that  the  level  of  the  water 
surface  in  the  channel  equals  the  elevation  of  the  channel  banks.  This  is 
obviously  not  true,  since  we  are  dealing  with  baseflow,  therefore  low  flows, 
for  the  purpose  of  groundwater  initialization.  This  is  a  point  where  more 
realistic  assumptions  are  needed  in  future  improvements  of  the  model. 
We  expect  that  given  the  approximation  used,  the  extent  of  the  near¬ 
channel  saturated  areas  are  overestimated  in  our  model. 
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3JE  FLOW  ROUTING 


i)  Overland  flow 


Runoff  routing  over  the  hillslope  considers  overland  flow  as  a 
quickly  channelized  process  rather  than  as  a  sheet  flow  regime  (Freeze, 
1980).  Therefore,  the  runoff  generated  in  an  element  is  not  allowed  to 
infiltrate  in  unsaturated  elements  located  downslope  from  it.  The  quality 
of  this  approximation  is  of  course  best  for  basins  of  high  drainage  density. 
Possible  infiltration  losses  in  microchannels  are  also  neglected.  Thus,  all 
runoff  generated  at  any  hillslope  element  is  routed  to  the  stream  network 
and  contributes  to  the  storm  hydrograph. 

Overland  flow  is  routed  using  a  constant  velocity  parameter, 
obtained  through  calibration.  This  is  done  lacking  parameter  estimates 
for  more  physically-based  routing,  and  for  simplicity  of  computation. 


ii)  Stream  flow 


Flow  routing  in  the  channel  network  also  uses  a  constant  travel 
velocity.  This  involves  two  approximations;  (1)  streamflow  velocity  is 
constant  along  the  channel  network  at  a  given  instant  in  time;  and  (2) 
streamflow  velocity  at  any  given  point  of  the  network  is  constant  at  all 
times  (independently  of  discharge).  While  the  first  approximation  is 
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supported  by  considerable  research  (as  we  will  see  below),  the  second 
approximation  is  reasonable  only  within  a  relatively  narrow  range  of 
discharge  values. 

Leopold  and  Maddock  (1953)  studied  the  relation  between  streamflow 
velocity,  V,  and  discharge,  Q.  They  performed  velocity  measurements  at 
different  locations  along  the  channel  network  at  different  times. 
Discharge  varies  with  location  at  a  given  time  (increasing  downstream) 
and  varies  with  time  at  a  given  location  (according  to  some  distribution  of 
discharge  values).  These  researchers  attempted  to  relate  velocity  to 
discharge.  Their  experimental  data  were  fitted  by  an  exponential  relation, 

V  =  C  •  Q«  (3.7) 

where  C  is  a  constant  and  a  is  a  parameter.  Equation  (3.7)  predicts  that,  at 
any  given  time,  as  discharge  increase  downstream  in  the  channel,  flow 
velocity  also  increases;  and  that,  at  aby  given  location,  flow  velocity  will 
vary  in  response  to  discharge  fluctuations.  The  fitted  value  of  a  for  a  fixed 
time  and  various  locations  is  not  the  same  as  for  a  fixed  location  at 
different  times. 

For  a  given  time  and  different  locations,  Leopold  and  Maddock  found 
an  average  value  of  a  =  0.1  for  semi-arid  regions  in  the  U.S.  This  value  of 
a  is  small  enough  that  we  may  expect  a  linear  approximation  to  the  V-Q 
relation  to  provide  a  reasonable  estimate  of  total  travel  time  in  the  channel 
network.  This  result  is  in  agreement  with  that  of  Pilgrim  (1977).  This 
researcher  conducted  tracer  studies,  measuring  total  travel  times  in 
various  stream  reaches.  The  tracer  studies  revealed  that  travel  times 
were  approximately  linearly  related  to  channel  length,  implying  that 


84 


estimation  of  travel  times  can  be  achieved  with  a  uniform  velocity 
parameter.  Beven  (1979)  also  has  shown  that  channel  kinematic-wave 
velocities  may  be  approximated  by  a  constant  value,  even  though  flow 
velocities  depend  on  discharge  rate.  The  assumption  of  uniform  velocities 
in  stream  channels  at  any  moment  in  time  has  been  made  by  Surkan 
(1974)  and  Rodriguez-Iturbe  and  Valdes  (1979),  in  their  Geomorphological 
Unit  Hydrograph  concept.  Kirkby  (1976)  has  shown  that  catchment 
response  becomes  progressively  more  linear  for  large  catchments,  as 
channel  travel  times  increasingly  dominate  the  hydrograph  form.  The 
approximation  is  best  for  medium  to  high  discharge  rates,  which  are  of 
interest  to  us. 

For  a  given  location  at  different  times,  flow  velocity  at  a  fixed 
location  will  vary  in  response  to  discharge  fluctuations.  Leopold  and 
Maddock  found  an  average  value  of  a  =  0.34  for  the  same  regions.  Thus, 
the  value  of  a  for  a  fixed  location  with  varying  time  is  considerably  larger 
than  that  for  fixed  time  at  various  locations.  One  possible  explanation  for 
this  fact  may  be  that,  at  a  given  time,  flow  increase  downstream  is 
generally  accompanied  by  a  broadening  of  the  channel  cross-section, 
which  allows  increased  discharge  for  an  approximately  constant  velocity. 
An  increase  in  discharge  in  time,  however,  is  accompanied  by  an  increase 
in  velocity  at  a  given  location.  If  a  large  range  of  discharges  in  time  is 
involved,  a  constant  velocity  value  may  not  be  expected  to  provide  good 
estimates  of  travel  time  in  the  channels. 
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3.F  PARAMETER  SENSITIVITY  AND  CALIBRATION  PROCEDURE 


This  section  explores  the  sensitivity  of  model  results  to  parameter 
values,  rainfall  rate  and  initial  moisture  condition.  Sub-Section  3.F.1 
explores  sensitivity  of  infiltration-excess  runoff  and  of  lateral  subsurface 
discharge  in  one  element.  Sub-Section  3.F.2  illustrates  sensitivity  of 
saturation-excess  runoff  to  terrain  morphology.  Sub-Section  3.F.3  deals 
with  calibration  of  conductivity  parameter  f. 


3.F.1  SENSITIVITY  IN  ONE  ELEMENT 


i)  Sensitivity  of  infiltration-excess  runoff  in  one  element  to  soil 
parameters  and  initial  moisture  condition 


In  this  section  we  investigate  model  sensitivity  to  soil  parameters 
Kq,  f,  9S,  0r  and  e,  rainfall  rate  R  and  initial  recharge  rate  Rj.  We  compare 

model  results  for  a  horizontal  hillslope  element  in  terms  of  ponding  time 
and  infiltration-excess  runoff  generated. 

Figures  3.1  through  3.5  illustrate  model  sensitivity  to  soil 
conductivity  parameters  K0  and  f.  We  see  that  larger  runoff  rates  result 
with  lower  values  of  K0  and  higher  values  of  f.  When  rainfall  rate  is  lower 
than  Kq,  runoff  rate  has  a  step  increase  from  zero  at  the  time  of  ponding, 
that  is,  when  Zt  reaches  zero.  Note  that  for  the  largest  values  of  f 
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considered,  f  =  0.01  mm'1,  and  f  =  0.005  mm*1,  runoff  is  generated  in  less 
than  24  hours  even  with  Kq  =  100  mm/hr,  a  value  ten  times  higher  than  the 

rainfall  rate. 

Figures  3.6  and  3.7  summarize  cumulative  runoff  and  ponding 

time,  respectively,  in  24  hours  under  a  constant  rainfall  rate  of  10  mm/hr, 

as  a  function  of  Kq  and  f.  Figure  3.6  shows  an  increase  in  cumulative 

runoff  with  f,  and  reveals  that  model  sensitivity  to  Kq  is  higher  for  lower 

values  of  f.  This  results  from  the  fact  that  if  soil  conductivity  decreases 

very  gradually  with  depth  (low  f  values),  soil  conductivity  at  the  wetting 

front  is  not  very  different  from  surface  conductivity,  and  the  rate  of 

infiltration  after  ponding  (i.e.,  after  we  have  Zt  =  0)  is  close  to  the  surface 

saturated  conductivity  Kq.  For  the  same  reason,  ponding  time  is  also  more 

sensitive  to  Kq  for  high  f  values,  as  shown  in  Figure  3.7. 

Figure  3.8  illustrates  model  sensitivity  to  rainfall  rate  R  for  a  soil 

% 

with  Kq  =  15  mm/hr.  The  difference  in  ponding  times  for  rainfall  rates 
lower  than  Kq  is  to  be  noted.  Figure  3.9  illustrates  model  sensitivity  to  the 
initial  recharge  rate  Rj.  Runoff  increases  with  Rj. 

Figures  3.10,  3.11  and  3.12  illustrate  model  sensitivity  to  porosity 
parameters  0S,  0r  and  e,  respectively.  These  parameters  are  used  in  the 

Brooks-Corey  parameterization  of  unsaturated  conductivity.  Predicted 
runoff  decreases  with  0S  and  increases  with  both  0r  and  e. 
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ETfiURE  3.1;  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  soil  conductivity  parameter  f.  Parameters  and 
variables  used  are  Kq  =  5  mm/hr,  0S  =  0.5,  0r  =  0.05,  £  =  4,  Rj  = 

0.1  mm/hr  and  R  =  10  mm/hr. 
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FIGURE  3.2;  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  soil  conductivity  parameter  f.  Parameters  and 
variables  used  are  Kq  =  10  mm/hr,  9S  =  0.5,  9r  =  0.05,  e  =  4,  R*  s 

0.1  mirvhr  and  R  =  10  mm/hr 
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figure  ZJh  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  soil  conductivity  parameter  f.  Parameters  and 
variables  used  are  Kq  =  15  mm/hr,  6S  =  0.5,  9r  =  0.05,  e  =  4,  R*  = 

0.1  mm/hr  and  R  =  10  mm/hr. 
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HGIIBEJLi.  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  soil  conductivity  parameter  f.  Parameters  and 
variables  used  are  Kq  =  20  mm/hr,  es  =  0.5,  9r  =  0.05,  e  =  4,  R,  = 

0.1  mm/hr  and  R  =  10  mm/hr. 
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FIGURE  3*5*  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  soil  conductivity  parameter  f.  Parameters  and 

variables  used  are  Kq  =  100  mm/hr,  0S  =  0.5,  9r  =  0.05,  e  =  4,  R,  = 
0.1  mm/hr  and  R  =  10  mm/hr. 
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solid  line:  fsl-lQ'^mm'^ 

long-dashed  line:  f=5-10'^mm'^ 

dotted  line:  f=l-10‘^mra  ' 1 

dot-dashed  line:  f=5-10"^mm'x 


FIGURE  3.R:  Sensitivity  of  the  cumulative  runoff  rate  after  24 
hours  to  conductivity  parameters  Kq  and  f.  Parameters  and 
variables  used  are  9S  =  0.5,  9r  =  0.05,  e  =  4,  R,  =  0.1  mm/hr,  and 


R  =  10  mm/hr. 
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FTGTTRE  3.7:  Sensitivity  of  the  ponding  time  to  conductivity 
parameters  Kq  and  f.  Parameters  and  variables  used  are  0S  = 
0.5,  0r  =  0.05,  e  =  4,  Rj  =  0.1  mm/hr,  and  R  =  10  mm/hr. 
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Runoff  rote  [mm/hr] 


Time  [hr] 


solid  line:  R  =  20  mnvhr 

long-dashed  line:  R  =  15  mm/hr 

dotted  line:  R  =  10  mm/hr 

dot-dashed  line:  R  s  5  mm/hr 

FIGURE  3.8:  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  rainfall  rate  R.  Parameters  and  variables  used  are  Kq 
=  15  mm/hr,  f  =  IQ^mnr1,  9S  =  0.5,  0r  =  0.05,  £  =  4,  and  R*  =  0.1 
mm/hr. 


95 


Runod 


solid  line: 
long-dashed  line: 
dotted  line: 
dot-dashed  line: 


Ri 

Ri 

Ri 

Ri 


=  1.0  mm/hr 
=  0.5  mm/hr 
=  0.1  mnvhr 
=  0.01  mm/hr 


EIGUBEJLSi  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  the  initiai  recharge  rate  Rj.  Parameters  and  variables 
used  are  Kq  =  15  mm/hr,  f  =  lO-amm*1,  9S  =  0.5,  9r  =  0.05,  £  =  4, 
and  R  s  10  mm/hr. 
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solid  line: 
long-dashed  line: 
dotted  line: 
dot-dashed  line: 
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9s  =  0.6  mm/hr 


FIGURE  2.IQ;  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  the  porosity  parameter  9S.  Parameters  and  variables 
used  are  Kq  =  15  mm/hr,  f  =  10-3mm-\  9r  =  0.05,  £  =  4,  and  R  = 
10  mm/hr. 
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FKiTIRE  StlLSensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  the  porosity  parameter  9r.  Parameters  and  variables 
used  are  Kq  =  15  mm/hr,  f  =  10-3mm*1,  9S  =  0.5,  e  =  4,  E,  =  0.1 
mm/hr,  and  R  =  10  mm/hr. 


96 


Runofl  rale  [mm/hr] 


solid  line:  e  =  mm/hr 
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dot-dashed  line:  e  =  mm/hr 


FTf?UKE  2J&  Sensitivity  of  the  infiltration-excess  runoff  rate  over 
time  to  soil  parameter  e.  Parameters  and  variables  used  are 
K<)  =  15  mm/hr,  f  =  lO^mm*1,  0S  =  0.5,  0r  =  0.05,  Rj  =  0.1 
mm/hr,  and  R  s  10  mm/hr. 
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ii)  Sensitivity  of  subsurface  discharge  from  one  element  to  soil 
parameters  and  terrain  inclination 


In  this  section  we  investigate  the  sensitivity  of  lateral  subsurface 
discharge  from  a  hillslope  element  to  soil  parameters  and  terrain  slope. 
We  will  consider  a  constant  rainfall  rate,  R  =  10  mm/hr,  and  various 
values  of  soil  parameters  K0  and  f.  Table  3.1  summarizes  the  runs 

performed  and  the  figure  where  the  corresponding  results  are  presented. 


Run 

Ko[mm/hr] 

f  [mm-l] 

1 

10 

5.10-3 

3.13 

2 

10 

1.10*3 

3.14 

3 

10 

1.10^ 

3.15 

4 

20 

5.10-3 

3.16 

5 

100 

5.10-3 

3.17 

TABLE  3.1:  Summary  of  runs  performed  to 
illustrate  the  sensitivity  of  subsurface  discharge 
rate  to  soil  conductivity  parameters  Kq  and  f. 

Figures  3.13  through  3.17  represent  the  evolution  in  time  of  (a) 
variables  Zf  and  Zt,  and  (b)  rate  of  subsurface  discharge.  Subsurface 

discharge  rate  is  expressed  as  volumetric  discharge  rate  per  unit  surface 
area  of  soil,  i.  e.  in  dimensions  [L/T].  In  all  runs,  the  highest  subsurface 
discharge  rate  was  obtained  with  a  =  45°,  which  results  from  the  fact  that 
factor  (cos(a)-sin(a))  in  Equation  (3.3)  for  Qout  has  its  maximum  value  at 
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Discharge  rate  (mm/hr) 


Time  [hr] 


in  time  of  (a)  variables  Zf  and  Zt,  and  (b) 
lateral  subsurface  discharge  rate  for  different  values  of  terrain 
inclination  angle  a.  Parameters  and  variables  considered  are 
Ko  =  10  mm/hr.  f  =  MO*3  mm*1,  9S  =  0.5,  9r  =  0.05,  £  =  4.0,  R,  = 
0.1  mm/hr,  and  R  =  10  mm/hr. 
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Discharge  rate  [mm/hr] 


Discharge  rate  [mm/hr] 


MW] 


fthure  3 .16-.E  volution  in  time  of  (a)  variables  Zf  and  Zt,  and  (b) 
lateral  subsurface  discharge  rate  for  different  values  of  terrain 
inclination  angle  a.  Parameters  and  variables  considered  are 
Kq  =  20  mm/hr,  f  =  5-10-3  mm*1,  9S  =  0.5,  9r  =  0.05,  e  =  4.0,  R,  = 
0.1  mm/hr,  and  R  =  10  mm/hr. 
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FIGURE  3.17: Evolution  in  time  of  (a)  variables  Zf  and  Zt,  and  (b) 
lateral  subsurface  discharge  rate  for  different  values  of  terrain 
inclination  angle  ct.  Parameters  and  variables  considered  are 
Kq  =  100  mm/hr,  f  =  5-10-3  mm-1,  9S  =  0.5,  9r  =  0.05,  e  =  4.0,  R,  = 
0.1  mm/hr,  and  R  =  10  mnvhr. 


a  =  45°.  The  split  between  the  curves  for  Zt  and  Zf  (which  occurs  at  Z*(R), 

given  in  Equation  (2.15)),  indicates  the  development  of  a  zone  of  perched 
saturation.  For  K0  £  R  (runs  1,  2  and  3),  Zt  equals  zero  from  the  time  of 
initiation  of  rainfall.  For  Kq  >  R  (runs  4  and  5),  perched  saturation 
develops  only  after  some  time  of  continued  rainfall.  Once  saturation  has 
developed  within  the  soil,  lateral  subsurface  discharge  is  initiated. 

Comparison  of  runs  1,  2  and  3  (Figures  3.13,  3.14  and  3.15)  show 
sensitivity  of  subsurface  discharge  rate  to  parameter  f.  The  relation 
between  subsurface  discharge  rate  and  parameter  f  is  nonlinear,  as  is 
apparent  from  Equation  (3.3)  for  Qout.  Note  that  in  (3.3)  Zf  and  Zt  are  also 

determined  by  f  (see  equations  of  evolution  for  these  variables,  summarized 
in  Section  3.B).  The  lateral  discharge  rate  is  initially  higher  for  f  =  510*3 
mm'1  than  for  f  =  110'3  mm*1  (Figures  3.13  and  3.14),  but  after 
approximately  20  hours  of  rainfall  (for  a  »  45°),  the  opposite  is  true. 
Subsurface  discharge  rate  for  f  =  1-10‘4  mm-1  is  much  lower  than  that  for 
110*3  mm-1  within  24  hours. 

Comparison  of  runs  1,  4  and  5  (Figures  3.13,  3.16  and  3.17)  shows 
sensitivity  of  lateral  subsurface  discharge  to  conductivity  parameter  Kq 
(note  the  different  vertical  scales  in  the  three  plots).  For  Kq  =  20  mm/hr 
and  Kq  =  100  mm/hr  (Kq  >  R),  subsurface  discharge  is  initiated  only  after 
perched  saturation  has  developed,  i.  e.  when  Zf  reaches  Z*(R).  Subsurface 
discharge  rate  increases  in  time  as  the  perched  saturation  zone  grows. 
The  curve  of  subsurface  discharge  rate  is  initially  concave  but  becomes 
convex,  with  the  inflection  point  occurring  at  the  ponding  time  (Zt  =  0). 

It  is  important  to  note  the  small  magnitude  of  the  values  of  lateral 
subsurface  discharge  rate  relative  to  the  rainfall  rate,  R  =  10  mm/hr.  The 
highest  subsurface  discharge  rate  in  24  hours  was  obtained  in  run  5  for  a 


106 


=  45°  (Figure  3.17)  and  this  is  only  approximately  0.01  mm/hr,  that  is  1,000 
times  smaller  than  the  rainfall  rate.  Therefore,  we  should  expect  lateral 
subsurface  flows  during  a  rainfall  event  to  have  very  little  influence  in 
determining  ponding  times  and  infiltration-excess  runoff  generation  on  a 
hillslope.  We  must  keep  in  mind,  however,  that  this  is  under  the 
assumption  of  isotropic  soils.  This  may  not  be  the  case  when  lateral 
hydraulic  conductivities  are  higher  than  vertical  ones  (see  Appendix  A). 


3.F.2  SENSITIVITY  OF  SATURATION-EXCESS  RUNOFF  TO 
HILLSLOPE  MORPHOLOGY 


The  saturation-excess  runoff  mechanism  is  determined  by  water- 
table  elevations;  where  the  water-table  is  at  the  soil  surface  (waterlogged 
areas),  direct  runoff  occurs.  The  procedure  used  to  estimate  the  position  of 
the  water-table  is  described  in  Section  3.D.  Here  we  investigate  sensitivity 
cf  estimated  water-table  elevations  to  conductivity  parameters  Kq  and  f. 
and  to  recharge  rate  Rj.  For  the  purpose  of  illustration,  we  will  use  the 
DEM  of  the  Sieve  basin  in  Italy,  and  will  assign  to  this  basin  different 
parameters  and  recharge  rate  values.  These  values  are  fictitious  and  are 
used  only  for  illustration  purposes.  In  Chapter  4,  the  model  is  applied  to 
the  Sieve  basin  using  estimated  real  parameters.  Table  3.2  summarizes 
the  illustration  runs  performed  and  the  results  obtained. 

We  see  that  the  percentage  of  basin  saturation  at  steady-state  is  the 
same  in  runs  1,  5  and  8;  in  runs  2  and  6;  and  in  runs  4  and  7.  Runs  5  has 
a  conductivity  10  times  smaller  than  run  1  but  an  f  parameter  that  is  10 
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times  smaller  as  well.  Therefore,  the  total  saturated  transmissivity  of  the 
soil  column,  which  is  approximated  by  Ko/f  (see  Section  3.D),  is  the  same 

in  the  two  runs.  The  result  of  same  transmissivity  is  a  similar  position  of 
the  steady-state  water-table.  Run  8  has  a  value  of  f  that  is  10  times  larger 
than  that  in  run  1.  However,  the  recharge  rate  considered  is  also  10  times 
larger,  and  the  result  is  a  similar  position  of  the  steady-state  water-table. 
The  similarities  of  water-tables  for  runs  2  and  6,  and  runs  4  and  7  has  the 
same  explanation.  We  conclude  that  the  extent  of  the  saturated  areas  is 
constant  for  constant  values  of  the  ratios  Ko/f,  Ko/Ri  and  Rj/f. 


Run 

K0  [mm/hr] 

f 

R{  [mm/hr] 

%  Saturation 

1 

100 

10-2 

0.01 

95.6 

2 

100 

10-3 

0.01 

40.8 

3 

100 

10"4 

0.01 

25.3 

4 

10 

10-2 

0.01 

99.3 

5 

10 

lO-3 

0.01 

95.6 

6 

10 

10^ 

0.01 

40.8 

7 

100 

10-2 

0.1 

99.3 

8 

100 

lO'3 

0.1 

95.6 

9 

100 

KH 

0.1 

40.8 

TABLE  3.2:  Summary  of  runs  performed  to  illustrate  sensitivity 
of  water-table  position  and  the  extent  of  saturated  areas  to 
hillslope  morphology. 


Figures  3.18,  3.19,  and  3.20  represent  water-table  depth  for  runs  1,  2 


and  3,  respectively. 


FTftTHtE  3.18:  Map  of  depths  (meters  )  below  the  surface  of  the  water-table 
in  steady-state  with  the  recharge  rate  Rj  =  0.01,  for  soil  parameters  Kq  = 
100  mm/hr  and  f  =  10*2  mm-V 


109 


FIGURE  3.19:  Map  of  depths  (meters;  below  the  surface  of  the  water-table 
in  steady-state  with  the  recharge  rate  Rj  =  0.01,  for  soil  parameters  Kq  = 

100  mm/hr  and  f  =  10*3  mm-1. 
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FlfiUREJl  Map  of  depths  (meters)  below  the  surface  of  the  water-table 
in  steady-state  with  the  recharge  rate  Rj  =  0.01,  for  soil  parameters  Kg  = 

100  mm/hr  and  f  =  10"4  mm*1. 


3.F.3  CALIBRATION  OF  CONDUCTIVITY  PARAMETER  f 


Parameter  f  is  the  rate  of  exponential  decrease  with  depth  of  the 
saturated  hydraulic  conductivity,  as  expressed  by  the  soil 
parameterization  equation  repeated  below  as  Equation  (3.8), 

K(z)  =  Ko-e  fz  (3  8) 

Beven  (1982)  summarized  the  results  of  fitting  Equation  (3.8)  to  a 
variety  of  soils  whose  conductivity  had  been  measured  over  depth  by 
different  researchers.  The  values  of  f  providing  the  best  fit  to  the  various 
soils  ranged  from  10*3  mm"1  to  10‘2  mm*1. 

In  the  usual  situation,  however,  such  studies  are  not  available, 
especially  if  the  catchment  in  question  is  large  and  includes  a  great  variety 
of  soil  types.  In  that  situation,  soil  parameterization  must  be  made  based 
on  assumptions.  In  absence  of  field  measurements,  we  assume  that  a 
single  value  of  /represents  all  soil  types  in  the  basin.  Figure  (4.7)  in 
Chapter  4  shows  simulated  basin  discharge  over  time  starting  at  full 
saturation,  for  various  values  of  /,  in  the  model  application  to  the  Arno 
River  catchment.  We  see  from  this  Figure  that  the  maximum  discharge 
from  the  basin  varies  approximately  linearly  with  /.  This  is  explained  by 
the  fact  that  the  transmissivity  of  a  soil  with  conductivity  given  by  Equation 
(3.8)  is  approximated  by  Ko/f  (see  Section  3.D),  and  therefore  discharge, 
Qout,  from  a  saturated  element  of  slope  tan(a)  is  approximated  by 
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Qout=  tan(a)-y 


(3.9) 


If  the  range  of  basin  baseflow  values  is  known  (as  usually  is,  or  can 
be  estimated  from  typical  regional  values),  then  an  important  physical 
basis  imposes  an  upper  limit  on  the  value  of  f.  The  value  of  f  must  be  such 
that  the  maximum  baseflow  obtained  by  the  groundwater  model  be  equal 
or  higher  than  the  maximum  baseflow  observed  in  the  basin.  However,  if 
we  have  reason  to  believe  that  the  soil  parameterization  given  by  Equation 
(2.12)  applies  only  to  a  limited  soil  depth  rather  than  to  the  entire  soil 
profile,  then  the  groundwater  model  cannot  be  used  to  predict  the  initial 
water-table  profile.  The  water-table  position,  and  the  extent  of  the  zone  of 
saturated  soil,  must  then  be  estimated  through  field  measurements.  In 
this  case,  parameter  f  can  be  obtained  by  calibration  with  observed  flows. 
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Chapter  4 


MODEL  RESULTS:  APPLICATION  TO  THE  SIEVE  RIVER  BASIN 


This  chapter  presents  the  application  of  the  rainfall-runoff  model  to 
the  Sieve  River  basin  in  the  Tuscany  Region,  Italy.  The  hydro-geophysical 
characteristics  of  the  Sieve  are  reviewed  in  Section  4.A.  Section  4.B 
describes  the  pertinent  data  available  for  the  basin  and  data  processing. 
Calibration  of  conductivity  parameter  f  and  routing  parameters  Vs  and  V0 

is  dealt  with  in  Sections  4.C  and  4.D.  Finally,  in  Section  4.E  we  present 
results  of  predicted  basin  response  for  observed  storms. 

The  Sieve  basin  has  an  area  of  approximately  840  Km2  and  is  one  of 
the  major  sub-basins  of  the  Arno  catchment  (8,000  Km2).  The  Amo  River 
crosses  the  cities  of  Florence  and  Pisa  and  has  produced  severe  flooding  of 
these  two  cities  in  the  past,  endangering  their  populations  and  raising 
both  national  and  international  concern  for  the  preservation  of  the  artistic 
and  historical  patrimony  of  Florence.  This  concern  has  led  to  the  creation 
of  the  Arno  Project,  which  aims  at  developing  a  monitoring  system  of 
hydraulic  variables  over  the  Arno  basin  which  will  increase  lead  time  of 
hydrologic  forecasting  for  the  Arno.  The  technologies  considered  for  the 
monitoring  system  include  telemetering  ground  sensor  networks,  a 
meteorological  radar,  airborne  and  satellite  sensors,  and  telemetering 
streamgages.  It  is  expected  that  such  technologies  be  implemented  in  the 
near  future.  The  new  distributed  information  that  will  thus  be  made 
available  should  be  fully  utilized  by  a  distributed  rainfall- runoff  model 
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such  as  the  one  presented  here.  As  we  will  see,  we  have  encountered  data 
limitations  in  our  model  application  to  the  Sieve  which  hopefully  will  be 
reduced  in  the  future. 


4*A  THE  SIEVE  RIVER  BASIN 


The  Sieve  basin  is  bordered  to  the  North-East  by  the  Appennine 
mountains,  with  a  point  of  maximum  elevation  at  Falterona  mountain 
(1657  meters  above  sea  level);  to  the  North  and  West  it  is  bordered  by  a 
series  of  mountain  ranges  of  peak  elevations  of  800  to  900  meters.  The 
average  elevation  of  the  basin  is  470  meters  above  sea  level.  The  basin  has 
moderate  to  strong  relief,  as  is  apparent  from  the  three-dimensional  view 
given  in  Figure  4.1.  The  confluence  of  the  Sieve  with  the  main  course  of 
the  Amo  is  located  19  Km  upstream  from  the  city  of  Florence.  The  basin 
has  approximately  63,000  inhabitants  (1971). 

Due  to  its  exposure  and  high  elevation,  the  Sieve  basin  is  subject  to 
intense  orographic  precipitation  events,  as  a  result  of  the  interception  of 
the  dominant  winds  coming  from  the  Mediterranean  Sea.  It  is  a 
characteristic  of  the  Mediterranean  climate  that  the  rain  season  coincides 
with  the  cold  season.  This  situation  is  particularly  propitious  to  the 
generation  of  high  runoff  volumes,  given  that  the  lower  evaporation  rates 
in  the  Winter  result  in  high  levels  of  soil  saturation  and  in  reduced 
infiltration. 

Monthly  precipitation  varies  widely  throughout  the  year,  with 
maximums  in  October-November  and  in  February,  and  minimum  flows  in 
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Aspect  of  the  topography  of  the  Sieve  basin, 
obtained  from  the  DEM. 
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July  and  October.  February  is  characterized  by  the  highest  amplitude  of 
observed  values  of  precipitation.  In  some  years,  February  will  be  among 
the  months  with  maximum  yearly  precipitation,  in  other  years  among 
those  with  minimum  precipitation. 

Table  4.1  shows  the  probability  distribution  (values  with  exceedence 
probability  of  10%,  50%,  and  90%)  of  monthly  minimum  flows,  in  cubic 
meters  per  second,  observed  at  Fornacina. 


4.B  DATA  AVAILABLE  FOR  THE  SIEVE  BASIN 


This  section  summarizes  the  data  available  for  the  Sieve  basin;  the 
characteristics  of  the  digital  elevation  map,  soils  data,  river  stage  data; 
and  rainfall  data. 


4.B.1  DIGITAL  ELEVATION  MAP 


The  DEM  available  for  the  Sieve  basin  has  a  square  grid,  400  meters 
to  the  side,  with  elevation  given  in  integer  meters,  with  North-South  and 
East-West  coordinates.  The  DEM  was  provided  by  the  Italian  Military 
Geographic  Institute  (IGMI),  and  was  obtained  through  the  IGMI 
standard  procedure  of  digitalization  of  contours  of  25  meters  interval  from 
a  topographic  map  at  the  scale  of  1:25,000,  utilizing  an  electronic  raster 
scanner  (A.  Carrara  et  al.,  1987).  The  IGMI  procedure  consists  of 
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month 

P=90% 

P=50% 

P=10% 

Maxim. 

Jan 

3.0 

7.5 

15.0 

22.6 

Feb 

3.0 

9.0 

17.0 

20.5 

Mar 

4.0 

7.0 

21.0 

24.0 

KHI 

2.5 

6.5 

14.5 

20.0 

3.0 

5.5 

9.5 

14.5 

Jun 

1.0 

3.5 

5.0 

10.0 

Jul 

0.5 

1.0 

2.5 

3.0 

Aug 

0.4 

0.9 

1.5 

2.0 

mm 

0.5 

0.9 

1.5 

2.0 

Oct 

0.5 

1.0 

4.0 

5.5 

Nov 

1.0 

4.0 

10.0 

18.5 

Dec 

2.0 

7.0 

15.5 

25.0 

TABLE  4.1  :  Distribution  of  monthly  minimum  flows  (in  m3/s) 
at  Fomacina.  P  represents  the  probability  of  exceedence.  The 
column  on  the  right  indicates  maximum  observed  values. 
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overlying  the  square  grid  on  the  topographic  map  and  identifying  the 
points  of  intersection  of  the  grid  lines  with  the  topographic  contours.  Lines 
of  maximum  slope  are  obtained  by  connecting  these  intersection  points, 
and  finally  an  elevation  value  is  assigned  to  each  grid  cell  by  interpolation 
between  the  maximum  slope  lines  (E.  Caporali,  1988).  This  procedure  of 
DEM  acquisit  .  is  subject  to  a  variety  of  errors,  associated  with  the 
approximations  inherent  to  the  establishment  of  the  lines  of  maximum 
slope,  rounding  errors,  and  operator  errors.  However,  a  limited  number 
of  validation  studies  have  been  made,  comparing  digitized  and  observed 
elevations,  indicating  a  maximum  error  of±  3  meters  (E.  Caporali,  1988). 

The  resolution  of  400  meters  seems  lower  than  desirable  for 
appropriate  representation  of  terrain  morphology,  given  the  marked  relief 
of  the  basin.  We  expect  that  at  this  resolution  terrain  slopes  computed 
from  the  DEM  underestimate  true  slopes.  Also,  given  the  grid  size, 
automatic  generation  of  the  channel  network  is  not  able  to  reproduce  the 
high  density  found  for  the  digitized  network. 

The  DEM  was  processed  to  obtain  estimated  distributed  terrain 
slopes  and  an  automatically-generated  channel  network.  The  algorithms 
used  for  both  purposes  were  developed  by  Tarboton  (see  Tarboton,  1989). 

Extraction  of  a  channel  network  from  the  DEM  requires  a  threshold 
contributing  area  The  most  appropriate  threshold  to  use  depends  on  the 
DEM  grid  size,  the  observed  drainage  density,  and  the  purpose  for  which 
the  extracted  network  will  be  used.  Figure  4.2  represents  the  map- 
digitized  stream  network  of  the  Sieve  basin.  It  is  apparent  that  the  detail 
and  high  density  of  this  network  cannot  be  reproduced  by  the  400  meter 
grid.  In  fact,  the  large  grid  spacing  of  the  DEM  represents  the  most 
important  limitation  to  network  density.  Carld  et  al.  (1986),  working  with 
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ftgijre  4.2:  Channel  network  obtained  from  digitalization  of  map  blue 
lines. 
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the  same  DEM  for  the  Sieve,  found  that  a  threshold  area  corresponding  to  8 
elements  (1.28  K2)  produced  a  network  which  compares  very  favorably  with 
that  represented  in  maps  at  the  scale  of  1:200,000.  We  have  adopted  the 
threshold  value  of  8  elements  based  on  Carla's  findings.  The 
corresponding  channel  network  is  represented  in  Figure  4.3.  This 
network  has  1084  stream  elements  (the  total  number  of  elements  in  Figure 
4.2  is  5252).  The  drainage  density  (average  length  of  channel  per  unit 
basin  area)  of  the  generated  network  is  0.52  Km/Km2.  The  drainage 
density  obtained  from  mapped  blue  lines  is  13.2  Km/Km2  (Carld  et  al., 
1986).  Again,  we  can  not  expect  to  reproduce  observed  drainage  density 
with  the  400  m  grid.  The  maximum  Strahler  order  is  IV.  Table  4.2 
represents  the  number  and  length  of  channels  by  Strahler  order. 


Strahler  order 

Number  of 

channels 

Length  of 

channels  [Km] 

I 

63 

253 

II 

14 

87 

III 

2 

38 

IV 

1 

17 

TABLE  4.2:  Number  and  length  of  channels  by  Strahler  order, 
obtained  from  the  400m  DEM  for  the  Sieve  basin  using  a 
threshold  contributing  area  of  1.28  Km2  (8  elements). 
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ftottre  4.3:  Channel  network  generated  from  the  400  meter  grid  DEM, 
using  a  threshold  area  equal  to  8  elements  C1.28  Km2). 
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4.B.2  SOILS  DATA 


No  field  measurements  of  soil  parameters  are  available  for  the 
Sieve.  However,  a  map  of  soil  types  exists  from  which  estimated  physical 
parameters  were  derived  based  on  measurements  conducted  for  similar 
soils  in  other  basins  in  the  region,  and  based  on  typical  parameter  values 
by  soil  type  provided  in  hydrology  literature.  The  map  of  soil  types  is 
represented  in  Figure  4.4  . 

Ranges  of  values  of  saturated  conductivity  of  the  topsoil  were  derived 
from  measurements  conducted  with  soils  of  assumed  similar  type  and 
characteristics  from  basins  in  the  same  region.  The  ranges  of  conductivity 
values  is  large,  varying  over  two  orders  of  magnitude.  We  will  use  these 
values,  described  as  topsoil  saturated  conductivity,  as  estimates  for  our 
conductivity  parameter  K 0.  Given  the  ranges  of  values  provided,  we  will 
take  the  median  value  from  each  range  to  represent  our  best  estimate  of  Kq 

for  a  given  soil  type.  The  range  of  conductivity  values  provided  and  the 
resulting  estimate  of  Kq  are  given  in  Table  4.3  . 

Estimated  values  of  total  porosity  were  also  obtained  from  field 
studies  conducted  on  similar  soils.  We  will  take  these  values  as  our  best 
estimates  for  parameter  0S  in  the  Brooks-Corey  parameterization.  For 
parameters  9r  (residual  porosity)  and  e  (pore-size  distribution  parameter) 
in  the  same  parameterization  equation,  we  will  adopt  values  indicated  in 
hydrology  literature  as  frequent  in  each  soil  textural  type  (Bear,  1972, 
Entekhabi,  1988).  The  parameter  values  used  are  summarized  in  Table 
4.3. 
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FIGURE  4.4:  Map  of  soil  types  of  the  Sieve  basin.  The  horizontal 
bars  represent  the  pattern  code  from  soil  type  1  (in  the  bottom) 
to  soil  type  17  (on  top).  The  length  of  each  bar  is  proportional  to 
the  fraction  of  the  basin  represented  by  the  respective  soil  type. 
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Soil 
Type  # 


Texture 

Class 


FS-FLA 


FL-FA 


FL-FLA 


FL 


F-FA 


FS 


FS-FA 


FS-FA 


FS-FA 


F-FS-FA-FL 


F-FS 


FS 


Detritic 


K  [mm/hr] 

Ko 

[mm/hr] 

es 

®r 

2-41 

21.5 

0.53 

0.02 

0.2-7 

3.6 

0.52 

0.04 

0.01  -  0.5 

0.25 

0.48 

0.09 

0.01  -  0.5 

0.25 

0.48 

0.09 

20-70 

45 

0.56 

0.11 

1-50 

25.5 

0.56 

0.09 

0.2-7 

3.6 

0.53 

0.11 

0.2  - 10 

5.1 

0.49 

0.11 

0.2  -  33 

16.6 

0.52 

0.06 

2.7-41 

21.8 

0.48 

0.04 

0.2  -  41 

20.6 

0.52 

0.07 

0.2  -  41 

20.6 

0.52 

0.07 

0.2  -  41 

20.6 

0.52 

0.07 

0.2-41 

20.6 

0.50 

0.07 

2.7-41 

21.8 

0.49 

0.03 

2.7-41 

21.8 

0.48 

0.041 

40 

40 

0.25 

0.02 

TABLE  4.3  :  Summary  of  soil  types  and  parameters  used  in  the 
model. 
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4.B.3  RIVER  STAGE  DATA 


The  Sieve  river  is  gaged  at  Fomacina,  a  location  near  the  confluence 
of  the  Sieve  with  the  main  course  of  the  Amo  river.  Stage  recording  is  not 
continuous  in  time;  it  is  initiated  manually  some  time  after  the  onset  of  a 
rainstorm  event  perceived  by  the  operator  to  be  important.  Therefore, 
existing  stage  records  are  of  short  duration  and  do  not  cover  inter- storm 
periods.  A  stage  curve  exists  for  the  Fornacina  gage  which  allowed 
conversion  of  stage  readings  to  flow  rates. 


4.B.4  RAINFALL  DATA 


Rainfall  measurements  at  intervals  of  20  minutes  are  available  for 
the  period  from  1968  to  1985,  for  four  raingages  located  inside  or  in  the 
vicinity  of  the  Sieve  basin.  The  Thiessen  polygons  corresponding  to  each  of 
these  raingages  are  represented  in  Figure  4.5  .'  We  see  that  one  gage, 
Borgo  S.  Lorenzo,  covers  approximately  75%  of  the  basin  area.  The 
available  rainfall  record  therefore  does  not  provide  the  spatial  resolution 
necessary  to  take  best  advantage  of  distributed  rainfall-runoff  modeling. 

The  rainfall  record  contains  several  gaps  which  reduce  the  number 
of  storms  for  which  we  have  both  stage  and  rainfall  data  to  be  used  for 
calibration  runs. 

Hourly  rainfall  records  with  higher  spatial  resolution  are  available 
for  3  rainfall  events;  of  November  of  1987,  December  of  1981,  and  February 
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FTflURE  4.5:  Thiessen  polygons  corresponding  to  the  4  raingages  that 
recorded  most  rainfall  events.  Only  one  raingage  is  located  inside 
the  Sieve  basin. 
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of  1968.  The  storm  of  November  of  1987  was  recorded  by  7  raingages.  The 
corresponding  Thiessen  polygons  are  represented  in  Figure  4.6.  The 
storms  of  October  of  1981  and  February  of  1968  were  recorded  by  6 
raingages;  the  same  as  the  storm  of  November  of  1987,  except  for  the  gage 
of  Vetta  Le  Croci. 


4.C  CALIBRATION  OF  CONDUCTIVITY  PARAMETER  f 


A  uniform  value  of  f  will  be  used  throughout  the  basin,  given  that  we 
do  not  have  information  on  soil  depth  for  the  different  soil  types.  Table  4.1 
above  gives  the  distribution  of  monthly  minimum  flow  values.  The  values 
with  an  exceedence  probability  of  50%  are  interpreted  as  representing 
basin  baseflow,  i.  e.  as  not  including  any  component  of  direct  runoff. 
These  baseflow  values  vary  from  0.9  to  9.0  m3/s. 

Figure  4.7  represents  the  discharge  curves  corresponding  to  various 
values  of  f,  using  the  estimated  values  of  surface  conductivity  and  porosity 
parameters  in  Table  4.3.  The  runs  represented  in  the  figure  were  initiated 
at  full  saturation  of  the  basin  (with  the  groundwater  coinciding  with  the 
terrain  surface),  and  the  initial  discharge  corresponds  to  the  maximum 
discharge  allowed  by  the  value  of  f  used.  We  see  that  the  values  of  f 
smaller  than  10'4  mm*1  do  not  provide  groundwater  discharges 
comparable  to  the  observed  50%  exceedence-probability  minimum  monthly 
flows  of  the  Winter  months.  Even  at  full  saturation  of  the  basin,  the 
groundwater  discharge  obtained  with  f  =  10*3  mm*1  was  approximately  ten 
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Ponte  2  Olmo 
Calenzano 
Borgo  S.  Lorenzo 
Vetta  Le  Croci 
V  i  1  lore 

Nave  di  Rosano 
Consuma 


FTfiTFRE  4.6:  Thiessen  polygons  corresponding  to  the  raingages 
that  recotded  the  event  of  November,  1987. 
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Groundwater  Discharge  [m^/s] 


12 


FTOURE  4.7:  Curves  of  predicted  groundwater  discharge  from 
the  Sieve  basin  for  various  values  of  parameter  f. 
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times  smaller  than  the  50 %  exceedence  probability  minimum  flow  for  the 
month  of  February. 

Only  the  value  f  =  10*4  mm-1  can  account  for  basin  baseflow  at  full 
saturation  of  the  basin  (at  time  zero).  Therefore,  if  we  assume  the  soils  to 
be  isotropic,  we  must  not  choose  a  value  of  f  smaller  than  10'4  mm*1.  For 
this  reason,  we  chose  the  value  f  *  10*4  mm*1,  a  very  small  value  but  which 
represents  a  maximum  established  by  the  observed  baseflow  values. 

It  should  be  noted  that  if  we  considered  the  soil  to  be  anisotropic, 
with  conductivity  in  the  direction  parallel  to  the  surface  higher  than  in  the 
direction  normal  to  the  surface  according  to  an  anisotropy  ratio  ar  (see 

Appendix  A),  then  the  same  discharge  at  full  saturation  would  be  obtained 
for  a  value  of  f  aT  times  larger  than  10*4  mm*1.  As  an  example,  in  Figure 
4.8  we  show  the  discharge  curve  for  f  =  10*3  mm*1  and  ar  =  10.  Discharge 

at  time  zero,  when  the  basin  is  fully  saturated,  is  the  same  for  f  =  lO^mm*1 
and  ar  =  10  as  for  f  =  10*4  mm*1  and  ar  =  0.  However,  the  time  rate  of  decay 

of  the  discharge  is  higher  for  the  anisotropic  soil. 

We  did  not  consider  anisotropy  for  the  soils  of  the  Sieve  given  that  no 
data  are  available  on  this  characteristic.  An  anisotropy  ratio  could  of 
course  be  obtained  through  calibration  of  model  results,  but  that 
calibration  would  be  interdependent  with  that  for  parameter  f.  The 
resulting  degrees  of  freedom  in  parameter  fitting  would  result  in  little 
physical  meaning  for  parameters  f  and  ar. 
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FIGURE  4. Si  Curve  of  predicted  groundwater  discharge  from  the 
Sieve  basin  for  f  =  10*3  mm'1  and  an  soil  anisotropy  ratio  of  ar=10. 
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4.D  CALIBRATION  OF  ROUTING  PARAMETERS  VSANDV0 


We  saw  in  Section  3.E  that  routing  parameters  Vs  and  V0  represent 
the  velocities  of  streamflow  and  overland  flow,  respectively,  assumed 
independent  of  flow  levels  and  spatially  uniform.  We  used  various  values 
of  these  two  parameters  and  in  different  combinations,  to  obtain  predicted 
hydrographs  for  the  calibration  storms.  The  best  combinations  of  Vs  and 
V0  values  for  each  storm  are  summarized  in  Table  4.4.  To  illustrate  the 

sensitivity  of  results  to  the  two  routing  parameters,  we  present  in  Figure 
4.9  the  predicted  hydrographs  obtained  for  different  values  of  Vs  and  V0  for 

the  storm  of  December,  1981. 

In  general,  the  sensitivity  of  results  to  routing  parameter  V0  was 
very  small  for  values  higher  than  0.4  Km/hr.  With  V0  lower  than 
0.4  Km/hr,  however,  hydrograph  dispersion  would  increase  markedly, 
with  the  shape  and  peak  intensity  not  being  preserved.  For  this  reason,  we 
considered  VQ  =  0.4  Km/hr  to  be  an  appropriate  value.  For  the  stream 
velocity,  V8  =  C.O  Km/hr  is  the  most  frequent  value  providing  the  best  fit. 


4B  RESULTS  FOR  OBSERVED  STORMS 


The  calibrated  parameter  values  that  we  will  use  are  Vs  =  6.0 
Km/hr,  V0  =  0.4  Km/hr,  and  f  =  10'4  mm*1.  In  this  section  we  present  the 
predicted  hydrographs  for  twelve  observed  rainfall  events  in  the  Sieve 
basin. 
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Storm 

Vs  [Km/hr] 

V0  [Km/hr] 

Feb  1968 

6.0 

0.4 

Dec  1968 

6.0 

0.4 

Jan  1969 

6.0 

0.4 

Dec  1975 

6.0 

0.2 

Dec  1976 

5.5 

0.4 

Feb  1977 

7.0 

0.4 

Jan  1979 

6.0 

0.3 

Dec  1981 

8.0 

0.4 

Nov  1982 

4.0 

0.4 

Feb  1983 

6.0 

0.4 

Jan  1985 

5.0 

0.4 

Nov  1987 

5.5 

0.4 

table  4.4;  Velocity  parameters  providing  the  best  fit  between 
predicted  and  observed  hydrographs  at  Fornacina  for  the 
rainfall  events  indicated. 


134 


0  to  20  30  *0  50  60  0  '0  20  30  *0  50  50 

Tim#  sinct  th#  initiation  of  romfall  {hr]  Tim*  smct  th*  initiation  or  ramtoil  [hrJ 


FIGURE  4.9  s  Sensitivity  of  the  predicted  hydrograph  to  the  velocity  parameters  V 
and  V0  for  the  storm  of  December,  1981.  Note  the  small  sensitivity  to  V0  for  V0  >  0. 
Km/hr.  The  figure  demonstrates  the  subjectivity  involved  in  selecting  one 
combination  of  Vs  and  V0  as  "the  best"  for  a  given  storm  (for  example,  not  the 
similarity  between  the  predicted  hydrographs  with  V9=  8  Km/hr  and  V0=  0.4  Km/hr 
and  with  Vs  =  7  Km/hr  and  V0  =  0.6  Km/hr). 
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Both  stage  and  rainfall  records  exist  for  these  twelve  rainstorms. 
All  events  originated  unusually  high  discharges.  Unfortunately,  no  data 
is  available  for  the  pre-storm  baseflow  for  any  of  the  events.  Therefore,  we 
cannot  determine  Rj  from  observed  baseflow  as  described  in  Section  3.D. 

In  face  of  the  unavailability  of  baseflow  values,  we  will  compute  predicted 
hydrographs  considering  three  hypothetical  pre-storm  baseflow  values, 
representing  a  typical  dry,  average  and  moist  condition  for  the  month  of 
the  event  in  question.  To  represent  the  dry  condition  we  will  use  the 
baseflow  value  that  has  an  exceedence  probability  of  90%  in  that  month;  for 
the  average  condition  we  will  use  the  value  with  50%  exceedence 
probability;  and  for  the  wet  condition  we  will  use  the  value  with  10% 
exceedence  probability.  As  examples,  Figures  4.10,  4.11  and  4.12  represent 
derived  steady-state  groundwater  positions  corresponding  to  the  wet, 
average,  and  dry  minimum  flow  for  the  month  of  November. 

These  results  are  presented  in  Figures  4.13  through  4.48.  By 
comparing  the  results  of  the  three  runs  corresponding  to  each  storm,  we 
see  that  for  some  storms  the  dry  initial  condition  provided  the  best 
hydrograph  prediction,  but  for  other  storms  the  average  or  the  wet  initial 
condition  was  best.  Using  the  minimum  deviation  between  predicted  and 
observed  peak-flows  as  the  criterion  to  define  the  best  good  prediction,  we 
see  from  the  Figures  that  5  storms  were  predicted  best  with  the  dry  initial 
groundwater  condition  (December  1968,  January  1969,  December  1975, 
November  1982,  November  1987),  4  with  the  average  condition  (February, 
1968,  February  1977,  January  1979,  January  1985  ),  and  3  with  the  wet 
condition  (December  1976,  December  1981,  February  1983).  Given  that  all 
12  storms  studied  had  high  observed  peak-flow  discharges,  it  is  somewhat 
surprising  that  in  5  of  them  the  dry  initialization  provided  the  best 
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estimated,  while  the  wet  initial  condition  was  best  in  only  3  cases.  A 
possible  explanation  for  this  is  the  fact,  already  noted  in  Chapter  3,  that  the 
model  to  estimate  the  depth  of  the  initial  water-table  tends  to  over-estimate 
the  extent  of  the  zones  of  saturation,  given  that  the  boundary  condition  for 
water-table  depth  adjacent  to  a  channel  is  approximated  by  the  elevation  of 
the  channel  banks.  If  the  extent  of  the  zones  of  saturation  are 
overestimated  for  a  given  water-table  recharge  rate,  then  a  better 
prediction  results  from  considering  a  lower  recharge  rate. 

To  investigate  a  relation  between  the  initial  condition  which  obtained 
the  best  estimate  and  the  actual  antecedent  condition  in  the  basin,  we  will 
use  an  antecedent  precipitation  index  as  an  indicator  of  the  degree  of 
wetness  of  the  basin.  Figure  4.49  represents  the  quality  of  the  prediction 
obtained  with  the  average  initialization  (50%  baseflow  exceedence 
probability),  in  terms  of  percentage  deviation  of  the  predicted  from  the 
observed  peak-flow,  against  an  antecedent  precipitation  index  (cumulative 
precipitation  over  the  preceding  30  days).  Only  9  of  the  12  storms  are 
represented  in  the  Figure  given  the  gaps  encountered  in  the  rainfall 
records  in  the  20-day  period  preceding  each  of  the  remaining  3  storms.  We 
should  expect  to  find  a  descending  tendency  in  the  plot,  i.  e.,  that  storms 
with  low  antecedent  precipitation  to  be  overpredicted  by  the  assumption  of 
an  average  initial  moisture  content  of  the  basin  and,  inversely,  storms 
with  high  antecedent  precipitation  to  be  underpredicted  by  that 
assumption.  This  expected  tendency  is  suggested  in  the  plot. 
Nevertheless,  the  small  number  of  points  in  the  plot  makes  it  difficult  to 
derive  a  definite  conclusion  about  any  tendency  in  the  results.  A  tendency 
would  mean  that  the  best  value  of  Rj  to  use  in  groundwater  initialization 
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during  operational  use  of  the  model  can  be  obtained  from  the  observed 
baseflow  at  the  storm  start,  as  we  suggest. 

We  consider  that  the  model  with  the  calibrated  f  and  velocity 
parameters  provided  overall  reasonable  estimates  of  peak  flow  rate, 
timeliness  of  peak  flow,  and  hydrograph  shape  in  general.  It  must  be 
noted,  however,  that  these  are  not  true  verification  runs,  as  the  storms 
utilized  here  are  the  same  used  for  calibration  of  the  routing  parameters. 
However,  the  velocity  parameters  providing  the  best  fit  for  the  different 
storms  showed  clear  consistency,  which  in  itself  provides  credibility. 
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FUTURE  4.10:  Depth  below  the  surface  (in  meters)  of  the  water-table  in  a 
dry  month  of  November.  The  water-table  is  in  steady-state  with  the 
basin  discharge  rate,  =  1.0  m3/s  (the  corresponding  recharge  rate  is 
Ri  =  Qfc /A  =  0.0043  mm/hr),  that  has  an  exceedence  probability  of  90%  in 
November. 


TTfiTTRF  4.11:  Depth  below  the  surface  (in  meters)  of  the  water-table  in  a 
average  month  of  November.  The  water-table  is  in  steady-state  with  the 
basin  discharge  rate,  Qb  =  4.0  m3/s  (the  corresponding  recharge  rate  is 
Rj  =  Qb/A  =  0.0171  mm/hr),  that  has  an  exceedence  probability  of  50%  in 

November. 
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FIGURE  4.12;  Depth  below  the  surface  (in  meters)  of  the  water-table  in  a 
wet  month  of  November.  The  water-table  is  in  steady-state  with  the 
basin  discharge  rate,  =  10.0  m^/s  (the  corresponding  recharge  rate  is 
Ri  =  Qb/A  =  0.0428  mm/hr),  that  has  an  exceedence  probability  of  10%  in 
November. 
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Storm  of  February,  1968 


Pre-storm  oaser'low  has  an  exceedence  probability  of  90% 


Time  since  the  initiation  of  ra infail  [hr] 


PURE  4.24:  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1968.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  90%)  in  the  month  of  February, 
Rj  as  0.0129  mm/hr. 
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Storm  of  February,  1968 


?re— storm  oasef'.ow  nas  an  exceeaence  prcDODtlity  of  5C% 

600  -r - - - - - 


0  10  20  30  4-0  50  60  70  30  90 

Time  since  the  initiation  of  rainfall  [hr] 


FIGURE  4.14:  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1968.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 

has  an  exceedence  probability  of  50%)  in  the  month  of  February, 
Rj  =  0.0386  mm/hr. 
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Flow  [m^/sj 


Storm  of  February,  1968 


Time  s:nce  :he  initiation  of  rcinfail  [hr] 


FIGURE  4.15;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1968.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  10%)  in  the  month  of  February, 
Rj  =  0.0728  mm/hr. 


Storm  of  December,  1968 


FIGURE 4.1  fi;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1968.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  90%)  in  the  month  of  December, 
Rj  =  0.0087  mm/hr. 


Storm  of  Decemoer,  1968 


Pre  —  scarm  oasefiow  nas  an  exceedence  proDabiiity  of  5C?o 


Mme  since  the  initiation  of  rainfall  [hr] 


FIGURE 4.17;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1968.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 

has  an  exceedence  probability  of  50%)  in  the  month  of  December, 
Rj  =  0.0300  mm/hr. 
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FIGURE 4.18:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1968.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  10%)  in  the  month  of  December, 
Ri  =  0.0664  mm/hr. 
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Flow 


Storm  of  January,  1969 


Pre-storm  oaser'low  has  an  exceedence  proDQDility  of  90% 
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FIGURE 4.19:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  January,  1969.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  90%)  in  the  month  of  January, 
Rj  =  0.0129  mm/hr. 
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Storm  of  January,  1969 


Pre— storm  oasefiow  nas  an  exceedence  oroDcoiiitv  of  50% 
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FIGURE 4.2Q;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1969.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 

has  an  exceedence  probability  of  50%)  in  the  month  of  January, 
Rj  =  0.0321  mm/hr. 
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Storm  of  January,  1969 


Pre— storm  casefiow  nas  an  exceedence  probability  of  10% 
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ftoitke  4.21 :  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1969.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  10%)  in  the  month  of  January, 
R|  =  0.0643  mm/hr. 
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Storm  of  December,  1975 


Pre— storm  Daseflow  nas  an  exceedence  probaDility  of  90% 


Time  since  the  initiation  of  rainfcli  [hr] 


FIGURE 4.22;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1975.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  90%)  in  the  month  of  December, 
Rj  =  0.0087  mm/hr. 


151 


Storm  of  December,  1975 


Pre— scorm  ocser'iow  nas  an  exceedence  proDcdility  of  50% 


Time  since  the  initiation  of  rainfall  [hrj 


FIGURE 4.23:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1975.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 
has  an  exceedence  probability  of  50%)  in  the  month  of  December, 
R|  =  0.0300  mm/hr. 
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iime  since  the  initiation  or  re inTCii  l_n rj 


fhturk 4.24:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1975.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  10%)  in  the  month  of  December, 
Ri  s  0.0664  mm/hr. 


Flow  [m^/s] 


Storm  of  December,  1976 


pre_storm  oaser'low  has  an  exceedence  prodaaility  of  90% 


FTGURE 4.25:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1976.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  90%)  in  the  month  of  December, 
R|  =  0.0087  mm/hr. 
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FIGURE 4.2ft.  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1976.  Groundwater  was 

initialized  m  steady-state  with  the  average  recharge  rate  (i.  e.  that 

has  an  exceedence  probability  of  50%)  in  the  month  of  December, 
Ri  =  0.0300  mm/hr. 
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Flow  [m^/bl 


Storm  of  December,  1976 


Pre— storm  oasefiow  has  an  exceedence  probability  of  \Q% 


FIGURE 4.27:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1976.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  10%)  in  the  month  of  December, 
Rj  =  0.0664  mm/hr. 
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Storm  or'  Febru.arq,  19  7/ 

J 

Pre-storm  ocsefiow  nas  an  exceed-:.' ^  probability  of  30% 
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Time  since  the  initiation  of  rainfall  [hr] 


EKfllRE 4,28;  Observed  (solid  line)  and  predicted  (dashed  hnej 
hydrographs  for  the  storm  of  February,  1977.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  90%)  in  the  month  of  February, 
Rj  =  0.0129  mra/hr. 
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Storm  of  February,  1977 


FIGURE 4.29:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  February,  1977.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 
has  an  exceedence  probability  of  50 %)  in  the  month  of  February, 
Rj  =  0.0386  mm/hr. 


158 


Storm  of  February,  197  7 


Pre-storm  Dasefiow  nas  an  exceeding  probability  of  10% 
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EK3UBJE 4.3Q;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1977.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  10%)  in  the  month  of  February, 
Rj  =  0.00428  mm/hr. 
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Storm  cf  January,  1979 


-'me  since  initiation  of  rcinfan 


FIGURE  4.31;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1979.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  90%)  in  the  month  of  January, 
Rj  =  0.0129  mm/hr. 


160 


Flow  [m-V b  | 


Storm  of  January,  1979 


■^re  —  sicrn  oasefiow  nas  an  exceedence  proDaoiiity  cf  50% 


Hme  sines  :he  initiation  of  reinfeii  F-ri 


FIGURE 1.32;  Observed  (solid  line)  and  predicted  (dashed  line 
hydrographs  for  the  storm  of  January,  1979.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 

has  an  exceedence  probability  of  50%)  in  the  month  of  January 
R|  =  0.0321  mm/hr 
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Storm  of  January,  1979 


Pre— storm  Daseflow  nas  an  exceedence  probability  of  10% 


FIGURE  4.33;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1979.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  10%)  in  the  month  of  January, 
Rj  ^  0.0643  mm/hr. 
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Storm  of  Decemoer,  1981 


Pre— scorn  ocsenow  nas  an  exceedence  oroDaDiiitv  of  90% 

600  -i - — - — - 1 - 


FIGURE 4.34;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1981.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  90%)  in  the  month  of  December, 
Rj  =  0.0087  mm/hr. 
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Storm  of  December,  1981 


pre  —  storm  ocsefiow  nas  an  exceedence  praoaaiiity  of  5C?o 
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FT  (PURE  4.35:  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  December,  1981.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 
has  an  exceedence  probability  of  50%)  in  the  month  of  December, 
Rj  =  0.0300  mm/hr. 
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Storm  of  December,  1981 
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Pre— storm  Daseflow  has  an  exceeaence  probability  of  10% 


Time  since  the  initiation  of  rainfall  [hr] 


FIGURE  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  December,  1981.  Groundwater  wag 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  10%)  in  the  month  of  December, 
Rj  =  0.0664  mm/hr 
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Storm  of  November,  1982 


Pre  storm  oasenow  nas  an  exceeaence  probaDiiity  of  90% 
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FIGURE 4.37;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  November,  1982.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  90%)  in  the  month  of  November, 
Rj  =  0.00428  mm/hr. 
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Storm  of  November,  1982 


Storm  of  February,  1983 


?re— storm  oaseflow  nas  cr«  exceeaence  probaDility  of  90% 
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FIGURE 4.40;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1983.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  90%)  in  the  month  of  February, 
Rj  =  0.0129  mm/hr. 
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Storm  of  February,  1983 


initial  recr.arge  rate  nas  a  5C%  exceedence  crooability 
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ftoitre 4.41;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1983.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 
has  an  exceedence  probability  of  50%)  in  the  month  of  February, 
=  0.0386  mm/hr. 


Storm  of  February,  1983 


Pre— storm  oaseflow  nas  an  exceedence  probability  of  10% 


Time  since  ;r.e  initiation  of  rainfall  [hr] 


FIGURE 4.42:  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  February,  1983.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  10%)  in  the  month  of  February, 
Rj  =  0.00428  mm/hr. 
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Storm  of  January,  1985 


Pre— scorm  oasefiow  nas  cn  exceedence  probability  of  90% 


Time  since  the  initiation  of  rainfcil  fhrl 


FIGURE 4.43:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  January,  1985.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  90%)  in  the  month  of  January, 
Rj  =  0.0129  mm/hr. 
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Storm  of  January,  ’985 


Time  since  me  initiation  of  rcinr'ctl  i  nr: 


FIGURE 4.44;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1985.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 
has  an  exceedence  probability  of  50%)  in  the  month  of  January, 
R,  =  0.0321  mm/hr. 
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Storm  of  January,  1985 


Pre— storm  oasefiow  nas  an  exceedence  proDcciiity  of  10% 
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FIGURE.4.45;  Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  January,  1985.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 

exceedence  probability  of  10%)  in  the  month  of  January, 
Rj  =  0.0643  mm/hr. 
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Storm  of  NovemDer,  1987 


Pre-srarm  oaseflow  nas  an  exceeaencs  proDGDiiity  of  2C% 
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FIGURE 4.46:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  November,  1987.  Groundwater  was 
initialized  in  steady-state  with  the  dry  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  90%)  in  the  month  of  November, 
Rj  =  0.00428  mm/hr. 


Storm  cf  Novemoer,  1987 


EIG.URE Observed  (solid  line)  and  predicted  (dashed  line) 
hydrographs  for  the  storm  of  November,  1987.  Groundwater  was 
initialized  in  steady-state  with  the  average  recharge  rate  (i.  e.  that 
has  an  exceedence  probability  of  50%)  in  the  month  of  November, 
Rj  =  0.0171  mm/hr. 
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Storm  of  Novemoer,  1987 


Prg  — storm  oasefiow  has  an  exceedence  probability  of  1  C% 
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FIGURE 4.48:  Observed  (solid  line)  and  predicted  (dashed  line) 

hydrographs  for  the  storm  of  November,  1987.  Groundwater  was 
initialized  in  steady-state  with  the  wet  recharge  rate  (i.  e.  that  has  an 
exceedence  probability  of  10%)  in  the  month  of  November, 
R}  =  0.0428  mm/hr. 
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FTfflfRE  4i4i?i  Deviation  of  predicted  from  observed  peak-flow  versus 
cumulative  precipitation  in  the  30  days  preceding  the  storm.  An 
inverse  tendency  is  suggested  in  the  plot,  but  the  small  number  of 
points  makes  it  dificult  to  derive  a  definite  conclusion. 
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Chapter  5 


SUMMARY  AND  CONCLUSIONS 


5 A  SUMMARY  OF  RESULTS 


We  have  presented  a  distributed,  physically-based  model  of  runoff 
generation  in  a  catchment,  for  operational  use  in  flood  forecasting.  The 
model  accounts  for  both  the  infiltration-excess  and  saturation-excess 
mechanisms  of  runoff  production,  and  for  lateral  subsurface  flows.  The 
kinematic  approximation  is  used  to  model  infiltration  and  subsurface 
flows.  The  effect  of  local  terrain  slope  (obtained  from  the  DEM  of  the 
catchment)  on  subsurface  flows  and  the  development  of  areas  of  saturated 
soils  is  accounted  for.  The  model  uses  spatial  discretization  into 
rectangular  elements  which  correspond  to  the  DEM  grid. 

The  model  was  applied  to  the  Sieve  catchment,  in  Italy,  and  used  to 
predict  hydrographs  for  12  recorded  rainfall  events.  For  each  event  three 
different  initial  groundwater  conditions  were  considered,  corresponding  to 
a  typical  dry,  average,  and  wet  state  of  the  basin  in  the  month  in  question. 
In  some  of  the  storms,  the  dry  initial  groundwater  condition  provided  the 
best  prediction,  in  other  storms  the  average  initial  condition  performed 
best,  and  yet  in  others  the  wet  initial  condition  did  best.  A  dependence  of 
the  best  initial  condition  on  the  antecedent  precipitation  recorded  in  the 
preceding  30  days  was  investigated.  Given  the  small  number  of  storms, 
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the  investigation  was  not  conclusive.  However,  the  data  suggest  that  the 
hydrographs  of  storms  with  high  antecedent  precipitation  tended  to  be 
under-predicted  by  the  average  groundwater  initialization  (50% 
exceedence  probability),  which  storms  with  low  antecedent  precipitation 
tended  to  be  over-predicted  by  the  average  initialization. 

Given  the  data  limitations,  in  particular  regarding  distributed 
rainfall  data  (approximately  75%  of  the  basin  area  being  covered  by  a 
single  gage  in  most  of  the  storms  studied),  and  unavailability  of 
hydrograph  recession  curves  for  calibration  of  parameter  f,  the  predictions 
provided  by  the  model  appear  to  be  quite  satisfactory.  In  general,  the 
observed  hydrograph  was  comprehensed  between  the  hydrographs 
predicted  with  the  dry  and  wet  (90%  and  10%  exceedence  probabilities) 
initial  groundwater  conditions,  or  did  not  distant  themselves  far  from 
these  predictions. 

5.B  LIMITATIONS  OF  THE  MODEL 

The  model  presented  contains  limitations  of  two  different  types; 
those  associated  with  the  basic  approximations  and  assumptions  made, 
and  those  that  result  from  unavailability  of  data  required  to  calibrate 
and/or  initialize  the  model.  Both  types  of  limitations  will  be  summarized. 

Important  limitations  are  associated  with  the  intialization  of  the 
moisture  profile  above  the  water-table.  We  recall  that  the  initial 
unsaturated  moisture  profile  above  the  water-table  is  approximated  by  a 
function  0(Rj,z)  which  assumes  that  there  is  a  constant  rate  of 


percolation, Rj,  at  all  depths  z  above  the  water-table,  and  that  capillary 
forces  play  no  role  in  inter-storm  subsurface  unsaturated  flows  or  in 
shaping  the  unsaturated  moisture  profile.  Thus,  this  moisture 
intialization  does  not  account  for  the  effects  of  evapotranspiration  during 
inter-storm  periods,  nor  does  it  represent  the  profile  of  the  capillary  fringe. 
The  importance  of  the  capillary  fringe  in  runoff  generation  has  been 
addressed  by  various  researchers.  Abdul  et  al.  (1989)  have  demonstrated 
how  the  proximity  of  the  water-table  to  the  surface  results  in  extremely 
small  saturation-deficits,  with  a  very  fast  rise  of  water-table  level  occuring 
upon  a  very  small  amount  of  infiltration,  generating  both  saturation- 
excess  and  fast  groundwater  discharge. 

Sivapalan,  M.,  et  al.  (1988)  have  suggested  an  approximation  to  the 
initial  moisture  profile  that  assumes  a  situation  of  equilibrium  with  the 
water-table  position  (no  vertical  flow).  The  resulting  function  of  soil 
moisture  with  depth,  0(Zwt,z),  has  a  simple  form  and  may  be  incorporated 

into  the  model  presented  here  in  a  most  straightforward  manner.  The 
function  suggested  by  Sivapalan  et  al.  was  not  introduced  in  our  model 
given  that  it  requires  the  relation  between  matrix  potential,  y,  and 
moisture  content,  0,  and  this  is  general  not  available.  However,  should  the 
relation  y(0)  be  known,  incorporation  of  the  function  0(Zwt,z)  used  by 

Sivapalan  et  al.  would  be  desirable. 

We  have  treated  soil  porosity,  9S,  and  the  residual  moisture  content, 
9r,  to  be  constant  with  depth.  In  general,  however,  both  these  parameters 

may  be  expected  to  vary  with  depth,  and  to  be  correlated  with  hydraulic 
conductivity.  Beveu,  et  al.  (1980)  fitted  analytic  forms  of  0s(z)  and  0r(z)  to 

measured  parameter  values,  and  found  reasonable  fitting  functions.  If 
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such  functions  are  known,  the  model  may  benefit  from  their 
incorporation. 

Finally,  water-table  initialization  assumes  water-surface  elevations 
in  the  stream  channels  to  equal  the  elevation  of  the  banks.  This  boundary 
condition  will  generally  lead  to  over-prediction  of  water-table  levels  and  the 
extent  of  the  zones  of  saturated  soils  (i.  e.  where  the  water-table  is  at  the 
soil  surface). 


5.B  SUGGESTIONS  FOR  FUTURE  RESEARCH 


Some  aspects  of  the  model  which  have  been  touched  upon  but  not 
explored  in  detail  deserve  further  research.  The  most  important  of  these 
aspects  is  perhaps  soil  anisotropy.  In  Appendix  A  the  infiltration  model 
presented  in  Section  2.B  was  extended  to  an  anisotropic  soil  where 
hydraulic  conductivity  is  higher  in  the  parallel  than  in  the  normal 
direction.  Incorporation  of  the  anisotropy  ratio  into  the  computation  of 
water-table  position  involves  only  substitution  of  the  lateral  saturated 
conductivity,  Kop,  for  Ko.  The  effect  of  considering  an  anisotropy  ratio  was 

illustrated  in  Figure  4.8.  However,  further  studies  of  sensitivity  are 
required. 

In  the  brief  attention  dedicated  to  the  subject,  it  was  shown  that  the 
ratio  of  anisotropy  affects  the  rate  of  decay  of  the  predicted  recession  curve. 
Inversely,  we  suggest  that  analysis  of  observed  hydrograph  recession  be 
used  to  estimate  the  anisotropy  ratio,  ar.  Consideration  should  also  be 
given  to  the  possibility  that  ar  be  a  function  of  soil  depth. 
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Another  aspect  of  the  model  which  requires  improvement  is  the 
approximation  made  of  the  water-surface  elevation  in  the  stream  channels 
by  the  elevation  of  the  channel  banks.  While  this  approximation  may  be 
reasonable  on  a  flat  floodplain,  it  may  be  considerably  in  error  in  most 
instances.  If  information  is  not  available  on  water-surface  elevations, 
improved  approximations  are  desirable  that  relate  this  variable,  for 
example,  to  flow  contributing  area,  link  magnitude  or  Strahler  order. 

Finally,  water-table  elevations  estimated  by  the  model  should  be 
validated  against  observed  water-table  position.  Ideally,  field 
measurements  should  provide  data  on  water-table  depths  at  different 
times  of  the  year,  and  such  measurements  be  used  in  calibration  of 
parameters  f  and  ar. 
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APPENDIX  A;  KINEMATIC  INFILTRATION  IN  A  VERTICALLY 
HETEROGENEOUS,  ANISOTROPIC,  SLOPED  SOIL 


In  this  Appendix  we  expand  the  analysis  performed  in  Section  2.B  to 
an  anisotropic  soil.  Many  natural  soils  are  anisotropic  with  conductivity 
in  the  direction  normal  to  the  surface  less  than  in  the  parallel  direction,  a 
condition  termed  lateral  isotropy.  Here  we  consider  infiltration  in  a  soil 
where  saturated  hydraulic  conductivity  is  higher  in  the  p  direction  than  in 
the  n  direction  (see  Figure-caption  2.2  in  Section  2.B  for  definition  of 
coordinate  directions).  The  conductivity  in  any  direction  decreases 
exponentially  with  depth  in  the  vertical  direction,  z.  This  condition  is 
particularly  prone  to  the  generation  of  lateral  subsurface  flow,  given  thaf  it 
combines  two  features,  heterogeneity  and  anisotropy,  capable  of  laterally 
diverting  infiltration. 


i)  Soil  parameterization  and  moisture  initialization 


Saturated  hydraulic  conductivity  in  both  the  parallel  and  normal 
directions  decreases  exponentially  with  depth, 


Ks  (z)  =  Kq  •  e-f  z 

P  P 


(A.l) 


K,  (z)  =  Ko  •  e"f  E 

n  n 


(A.2) 
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where  KSp(z)  and  KSn(z)  are  the  saturated  conductivities  in  the  p  and 
n  directions  at  depth  z;  K0p  and  Kon  are  the  saturated  hydraulic 
conductivities  in  the  p  and  n  directions  at  the  soil  surface;  and  f  is  a 
parameter  of  dimension  Lr1. 


The  saturated  hydraulic  conductivities  in  the  p  and  n  directions  are 
related  through  the  dimensionless  anisotropy  ratio,  ar,  defined  as 

Kj  (z)  K0 

Md  =  Kf=a'>1  (A.3) 

n  n 


Assuming  the  same  Brooks-Corey  parameterization  in  all  directions 
then  the  unsaturated  hydraulic  conductivities  are 


(A.4) 


(A.5) 


and  we  have,  for  unsaturated  soils  as  well, 

Kp(e,z)  _  Kop 

KnfO.z)  ~  =  3r  (A.6) 


As  in  Section  2.B,  the  initial  moisture  condition  is  determined  by  an 
antecedent  rate  of  recharge,  Rj. 
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ii)  Unsaturated  infiltration 


The  wetting  front  is  unsaturated  for  rainfall  rates  lower  than  the 
surface  saturated  conductivity,  and  for  wetting  front  depths  smaller  than 
the  critical  depth  Z*(R).  All  equations  derived  under  point  ii)  are  for  z  '  Z r 
with  Zf  <  Z*(R).  Letting  K(9,z)=R  and  0=0S  in  (A.5)  and  solving  for  z,  we 
obtain  the  expression  for  Z*(R), 

Z*(R)  =  j  •  ln^— (A.7) 


For  an  anisotropic  soil  sloped  at  an  angle  a,  we  can  write  the 
expressions  for  the  components  of  unsaturated  flow  in  the  main  directions 
of  anisotropy,  by  considering  the  components  of  the  gravitational  gradient 
vector  Jz  =  1.0  in  those  directions, 


Jp  =  sin  (a) ; 

qp  =  Jp  •  Kp(0,z) 

(A.8) 

Jn  a  cos  (a) ; 

=  Jn  ■  Kn(0,z) 

(A.9) 

Soil  anisotropy  determines  that  flow  be  at  an  angle,  a",  with  the 
vertical  direction.  Figure  A.l  represents  the  geometrical  situation  from 
which  we  can  write  (using  Equations  (A.6),  (A.8)  and  (A. 9)), 
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FIGURE  A.i:  Components  of  flow  in  the  z 
and  h,  and  p  and  n  directions.  Flow  is  in 
the  direction  indicated  by  s. 
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fq  'i 

(a+a")  =  tan-1  — - 

l  On  ; 


=  tan  1(ar-tan(a)) 


and  solving  for  a", 


a"  =  tan  1(ar-tan(a))  -  a, 


(A.  10) 


Angle  a"  is  larger  for  higher  anisotropy  ratios  and  higher  terrain 
slopes.  Note  that  for  ar  =  1  (A.  10)  gives  a"  equal  to  zero.  That  is,  for  an 
isotropic  soil,  the  above  equations  give  unsaturated  flow  in  the  vertical 
direction,  as  expected. 

Angle  a"  is  needed  for  computation  of  the  unsaturated  flow  rate. 
Designating  the  direction  of  flow  by  s,  we  can  derive  an  expression  for  qs 
from  continuity  considerations.  Let  us  consider  rainfall  at  a  rate  R  over  a 
unit  horizontal  surface  area,  and  its  infiltration  into  the  soil.  Figure  A.2 
represents  a  cross-sectional  area.  From  continuity,  the  flow  per  unit 
width  perpendicular  to  segment  a  (which  has  unit  length)  must  equal  the 
flow  per  unit  width  perpendicular  to  segment  6  i.e., 

R  •  1.0  =  qs  •  b  (A.ll) 


From  geometrical  considerations,  we  can  write 
k  _  cos(a+a") 

cos(a)  (A.12) 

Substituting  (A.12)  into  (A.ll)  and  solving  for  qs  we  obtain 
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FIGURE  A.2:  Unsaturated  infiltration  in  a  laterally 
isotropic  soil. 
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(A.  13) 


_  cos(g) 
<*s"~  cos(a+a") 


We  will  define 

-  _  cos(a) 

^  cos(a+a") 

(§  >  1.0  for  ar  >  1.0) 

and  thus  we  can  write 

qs  =  $-R 


(A.  14) 


(A.15) 


From  qs  we  can  finally  derive,  from  geometrical  considerations,  the 
expressions  for  the  vertical  and  horizontal  components  of  flow, 

qz  =  cos(a")q3 

=  cos(a"H-R  (A.  16) 

qh  =  sinCa")^ 

=  sin(a")-^R  (A.17) 

The  lateral  component  of  flow  increases  with  a"  and  thus  is  larger 
for  higher  slopes  and  higher  anisotropy  ratios. 

We  can  find  0(R,z)  by  noting  that 
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cjn  =  cos(a+a")-qs 


=  cos(a)-R  (A.18) 

Equating  (A.18)  to  (A. 9)  we  obtain  Kn(9,z)  =  R.  Letting  Kn(9,z)  =  R  in 
(A.5)  and  solving  for  9(R,z)  we  obtain 

9(R,z)  =  j e  •  (9s-9r)  •  efz  +  0r  (A.19) 

Note  that  the  soil  saturates  at  the  surface  (9=9S)  for  a  rainfall  rate 
equal  to  the  surface  saturated  conductivity  in  the  normal  direction,  Kon. 

From  derivation  similar  to  that  for  (A.  16),  we  have  below  the  wetting 
front  (i.e.,  for  z  >  Zf), 


cfeCRi)  =  cos(a")-^R,  (A  2Q) 

Through  a  derivation  analogous  to  that  for  the  equation  of  evolution 
of  Zf  in  Section  2.B,  point  ii),  and  upon  substitution  of  (A.16)  and  (A.20),  we 

obtain 


dZf  R-  Rj 

di T  "  0(R,Zf)  -  0(R,,Zf) 


(A.21) 
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iii)  Saturated  infiltration 


Saturated  infiltration  occurs  when  either  the  rainfall  rate  is  higher 
than  the  surface  saturated  conductivity;  or  the  wetting  front  has 
penetrated  beyond  the  critical  depth  Z*(R).  The  equations  derived  under 
point  iii)  are  for  Zt<z<Zf. 

For  the  reasons  seen  in  2.B  based  on  continuity  considerations,  the 
normal  component  of  flow  within  the  saturated  zone  is  limited  by  flow  at 
the  wetting  front,  and  therefore  qn(z)  is  a  constant  for  Zt  <  z  <  Zf.  Thus, 

qn(z)  =  qn(Zf)  (A.22) 

At  the  wetting  front,  we  have  Jz(Zf)  =  1  and 

Jn(Zf)  =  cos(a)  Jz(Zf)  =  cos(a)  ; 

qn(Zf)  =  Jn(Zf)  •  Ksn(Zf)  =  cos(a)-K0n  e"f  Z' 


(A.23) 

(A. 24) 


Substituting  (A. 24)  in  (A.22),  we  obtain 

qn(z)  =  cos(a)-K0n'e’fZf  (A.25) 

The  potential  gradient,  and  the  flow,  in  the  parallel  direction  are 
derived  as  follows, 
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_  ,  3<X>(p,n)  3  f  ,  x,  P(n)  1 

Jp(z)s— 55 — me(p'n)+_rJ 


...  dn  3  P(n)  .  .  . 

=  sin(a)  -  — •  -r - — -  =  sin(a) 

dp  on  y 


and 


Qp(z)  =  Jp(z)-Kg  (z)  =  sin(a)*K(^-e 


-fz 


As  in  part  ii)  and  using  (A.25)  and  (A. 27),  we  can  write 


tan(a+a'(z))  = 


^  =  tan(a)areHzr‘> 

Qn(z) 


and  solving  for  a"(z), 


a'(z)  =  tan  \  tan(a)aref  (Zf  z)  ]  -  a  , 


From  geometrical  considerations,  we  have 


,  .  sin(a'(z))  ,  . 
qh(z)“  sin(a+a'(z))  'qP(z) 


Substituting  (A.27)  in  (A.29)  we  obtain  for  q^Cz) 


.  .  sin(a)-sin(a’(z))  -fz 
qi>(z)  =  sinCa+aTz))  K°p e 


(A.26) 


(A.27) 


(A.28) 


(A.29) 


(A.30) 
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From  geometrical  considerations,  we  have 


,  .  cos(a'(z))  .  , 
Z  ~  cos(a+a'(z))  qn  Z 


(A.31) 


Substituting  (A.25)  in  (A.32)  we  obtain  for  q^z) 

,  v  cosCal-cosCa'Cz))  Tr  _f-z( 

^ 2  “  cos(a+a’(z))  ^ 6  (A.32) 


At  the  wetting  front  (i.e.  making  z=Zf),  Equation  (A.28)  equals 
Equation  (A.  10),  i.e.  a’(Zf)  equals  the  angle  of  unsaturated  flow,  a",  which 
does  not  depend  on  z.  Thus,  for  z  =  Zf  (A.32)  becomes 


qz(Zf)  =  cos(a")  ^Kon  e-f  Zf 


(A.33) 


A  derivation  analogous  to  that  for  dZf/dt  in  Section  2.B  leads  to  the 
following,  upon  substitution  of  (A.33)  and  (A. 20) 


dZf  KoV^-Ri 

dt  ""  es  -  9(Ri,Zf)  (A-34) 


A  derivation  analogous  to  that  for  dZt/dt  in  Section  2.B  leads  to  the 
following,  upon  substitution  of  (A.33)  and  (A.  16) 


dZt  Ko;e"fZf-R 

dt  “  0,—  0(Rj,Zt)  (A-35) 
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iv)  Summary  of  computations  in  one  element 


The  general  expression  for  the  volumetric  discharge  from  an 
element  to  the  element  located  downslope  from  it  is 


(A.36) 


Upon  substitution  of  Equations  (A.  17)  and  (A.30),  (A.36)  becomes 


Qoot=w 


r2* 

sin(oc")£ 

Jo 


„  r,  ,  .  P  sin(a)-sm(a  (z))  _f.z  , 

•R-Zt-dz  +  - r—. - rr-rr — Ko  e  dz 

1  Jz^  sin(a+a  (z)) 


(A.37) 


Upon  substitution  of  Equation  (A.28)  for  a'(z)  for  Zt  £  z  <  Zf  in 
Equation  (A.37)  we  derive  (see  derivation  in  Appendix  B,  part  ii)). 


Qcut=W|  sin(a"K-R-Zt  + 


Kop-cos(a)-sin(a)-e  fZf- 


l.(ef<ZrZ-) 


-1)  -  (Zf-Zt) 


(A.38) 


Below  we  summarize  the  computations  in  one  element.  In  the 
equations  below,  A  designates  element  area.  Estimation  of  Q,n  is  achieved 

through  the  element  coupling  scheme  which  is  explained  in  Section  3.B. 


dMt  dZr 

=  —-eiR^Zf)  +  Min(R,K0  )  + 
at  dt  n 


Qin  “  Qout 
A 
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Mu  =  Mt-0s.(Zf-Zt) 

_i_ 

0(R,z)  =  •  (0s-0r)  •  e7'z  +  0r 

Qout=  w'[  sin(a’’)^-R-Zt  + 

Kc^-cos(a)-sin(a),e"fZf-  l  (ef{zfz‘)-l)  -  (Zf-Zp 

a"  =  tan_1(ar-tan(a))  -  a 
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APPENDIX  B:  SOLUTIONS  OF  THE  INTEGRALS  IN  THE 
EXPRESSION  FOR  LATERAL  DISCHARGE  FROM  AN  ELEMENT 


In  this  Appendix  we  present  the  solution  of  the  integrals  in 
Equations  (3.3)  and  (A.37)  leading  to  Equations  (3.3)  and  (A.38), 
respectively. 

i)  Solution  of  the  integral  in  Equation  (3.3)  for  lateral  discharge 
from  an  isotropic  element 


Equation  (3.3)  is  repeated  here  as  (B.l), 


Q-=wf 


•7 

^  tan(aHan(a'(z))  Tr  _f  z 
^  tan(a)+tan(a'(z))  °G  '  Z 


(B.l) 


The  expression  for  the  tangent  of  the  sum  of  two  angles,  a  and  p, 
(see  e.g.  Simmons,  1985)  is 


tan(a+P)  = 


tan(a)  +  tan(p) 

1  -  tan(a)-tan(P) 


(B.2) 


From  (B.2),  we  can  write 


tan(a+a'(z))  = 


tan(a)  +  tan(a'(z)) 
1  -  tan(a)-tan(a'(z)) 


(B.3) 


From  (B.3),  the  ratio  in  (B.l)  can  be  written  in  the  following  form, 


206 


tan(a)-tan(a'(z)) _ 1 _ 1 

tan(a)+tan(a'(z))  ~  tan(a)+tan(a'(z))  tan(a+a'(z)) 


(B.4) 


and  Equation  (B.l)  becomes 


Q, 


out 


.wfr _ i 

Jr  L  tan(a)+t£ 


^  L  tan(a)+tan(a’(z))  tan(a+a’(z))  J 


Kn-e_fz-dz 


or 


Qout =  WKq- 


I 


^  tan(a)+tan(a'(z)) 


e~fzdz 


•  t  1  c~fzd: 
4  tan(a+a’(z)) 


(B.5) 


To  solve  the  first  integral  in  (B.5)  we  must  substitute  an  expression 
for  tan(a’(z)).  Substituting  Equation  (2.25)  for  oc’(z),  we  have 


tan(a'(z))  =  tan[  tan  1(ef' (Zrz)-tan(a))  -  1  ] 


(B.6) 


From  expression  (B.2)  for  the  tangent  of  the  sum  of  two  angles,  (B.6) 
becomes 


tan(a'(z))  = 


ef  (Z|—  z)-tan(a)  -  tan(a) 
1  +  ef  (Zf  “  z)-tan2(a) 


(B.7) 
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Upon  substitution  of  (B.7),  the  first  integral  in  (B.5)  becomes 


^  .  ef(Zf  z)-tan(a)-tan(a) 

tan(a)  + - - - 

l+ef(zrz,'tan2(a) 


-e"f'z-dz 


Upon  simplification  of  the  integrand  expression  and  final 


integration, 


f  ^  1+  ef(2rz)-tan2(a) 

^  tan(a)-efZf-(tan2(a)+l) 


(Zf-ZJ  +  tan^a)- Y-(ef(Z^Z‘)-l) 
tan(a)efZf(tan2(a)+l) 


(B.8) 


Solving  Equation  (2.25)  for  (a+a’(z)),  we  obtain 


tan(a+a'(z))  =  ef  (Zrz)-tan(a) 


To  solve  the 


we  must  substitute  the 


expression  for  tan(a+a'(z)).  Substituting  the  above  expression  for 
tan(a+a'(z))  in  (B.5),  the  integral  becomes 


Jz  _f'(Zf-z)  . 


•e~fz,dz 


A  en  f~ZJ,tan(a) 


(B.9) 
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Upon  simplification  of  the  integrand  expression  and  final 
integration, 


_1 _ 

ftan(a) 


•dz 


2f-Zt 
ef  Zf  tan(a) 


(B.10) 


Substituting  (B.8)  and  (B.10)  for  the  two  integrals  in  (B.5),  and  upon 
simplification,  we  obtain  the  final  expression,  which  was  presented  also  as 
Equation  (3.3) 

0^=  W-Ko-tan(a)— •e'f'Zff  l(ef(z^Zt)-l)  -  (Zf-ZJ 
tanla)+l  L f  -I 

and  given  that  we  have  tan2(a)+l  = - ^ — ,  the  above  equation  becomes 

cosla) 


Qaut=  W-Kocos(a)sin(a)e 


-f-z, 


I.(ef(zrztLi)  _  (Zr-Zt)  (B.ll) 


ii)  Solution  of  the  integral  in  Equation  (A*37)  for  lateral  discharge 
from  an  anisotropic  element 


Equation  (A.37)  is  repeated  here  as  (B.12) 


+ 


f^sinCcQ-sinCa'Cz)) 

J  ^  sin(a+a'(z)) 


•Ko  e-f  z  dz 

p 


(B.12) 
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The  ratio  in  the  second  integral  of  (B.12)  is  in  fact  equal  to  the  ratio 
in  the  integral  of  (B.l),  as  shown  next.  Upon  substitution  of  the  formula 
for  the  sine  of  the  sum  of  two  angles  and  dividing  numerator  and 
denominator  by  (cos(a)-cos(a')), 

sin(a)-sin(a'(z)) _ sin(a)-sin(a'(z)) _ 

sin(a+a'(z))  sin(a)-cos(a'(z))  +  cos(a)-sin(a'(z)) 

„  tan(a)-tan(a'(z)) 

tan(a)  +  tan(a'(z))  (B.13) 

Given  the  similarity  between  the  second  integral  in  (B.12)  and  the 
integral  in  (B.l),  their  solution  is  alike  and  we  obtain  (B.14),  which  was 
given  as  Equation  (A.38), 

Qout=  W  [  sin(a"K-R  Zt  + 


K^-cos(a)-sin(a)-e  rZf-  j-(ef  (ZrZt)-l)  -  (Z^Z,) 


(B.14) 


APPENDIX  C:  DERIVATION  OF  THE  LIMITS  OF  POTENTIAL 
GRADIENTS  AND  LATERAL  DISCHARGE 


In  this  Appendix  we  derive  the  limits  of  the  potential  gradients  in 
the  vertical  and  horizontal  directions,  Jz(z)  and  J^Cz),  and  in  the  direction 
of  flow,  Js(z).  We  also  derive  the  limit  for  the  lateral  discharge,  Qout. 

The  equations  for  a'(z)  and  Jz(z)  and  Jh(z)  were  given  in  Equations 

(2.25),  (2.32)  and  (2.30),  respectively.  Here  we  write  these  three  equations 
for  z=Zt,  i.  e.  for  the  top  of  the  zone  of  perched  saturation. 


a(ZJ  =  tan  1[ef'(Zf'Zt)-tan(a)]  -  a 


(C.l) 


T  (7  v _ tan(a) 

2  ^  ~  tan(a)  +  tan(a’(Zt)) 


(C.2) 


Jh(ZJ  =  (1-J2(z))tan(a) 


(C. 


Under  a  constant  rainfall  rate  R,  Zt  decreases  as  Zf  increases  with 
time,  and  Zt  will  eventually  reach  the  soil  surface.  As  the  difference  (Zf  Zt) 
increases,  angle  a’(Zt)  increases.  For  very  large  f  and  (Zf  Zt),  Equation 
(C.l)  approaches  (90°-a),  i.  e.  flow  is  nearly  parallel  to  the  terrain  surface. 

HZf-Zt)  — >  oo;  Ot^Zf)  — >  (90°  —  Ct)  /fi 


From  trigonometry,  we  have  that 
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tan(90°-  a)  = 


tan(a) 


(C.5) 


Using  (C.4)  and  (C.5),  (C.2)  becomes, 


tan(a) 


tan(a)  + 


tan(a) 


Multiplying  numerator  and  denominator  by  tan(a)  and  upon 
simplification  we  obtain, 


f-(Zf — Z()  — >  Jz(Zt)  — >  sin  (ot) 


(C.6) 


Substituting  (C.6)  in  (C.3)  and  upon  simplification  we  obtain, 


f^Zf-Zt)  J^Z,)  cos(a)-sin(a) 


(C.7) 


For  Js(Zt)  we  have,  from  geometrical  considerations, 


J8(Zt)  =  VJh(Z/+Jz(Z/ 


(C.8) 


Substituting  (C.6)  and  (C.7)  in  (C.8)  we  obtain, 


f^Zf-Z,;)  ->  oo;  ^(Z,)  — >  Vsin4(a)  +  cos2(a)-sin2(a) 


which  upon  simplification  becomes, 
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f(Zr-Zt)  - *  °°> 


JgCZt)  sin(a) 


(C.9) 


The  lateral  subsurface  discharge  (for  an  isotropic  soil)  was  given  in 
Equation  (58),  repeated  here  as  (C.10), 

Qout=  W-Ko-cos(a)-sin(a)-e“fZ4  -kef(Zrz‘)  -  1)  -  (Zj-ZJ 

L f  J  (C.10) 


Equation  (C.10)  can  be  written  as, 

Qout=  W  K0  cos(a)-sin(a)-  ~(e'fZ*  -e_fZf )  -  •e-fZ'-(ZrZt) 


(C.ll) 


For  large  Zf  and  small  Zt,  Equation  (C.ll)  tends  to, 

Ko 

Zf  — >  oo,  Zt  0  ;  Qout^  W  —  cos(a)  sin(a) 


(C.12) 


We  saw  in  Section  III.D  that  total  transmissivity  is  approximated  by 
Ko/f  for  large  soil  depths.  Equation  (C.7)  gives  the  limit  for  the  lateral 
hydraulic  gradient  to  be  (cos(a)-sin(oc)).  Therefore,  Equation  (C.12)  equals 
the  limit  of  the  product  of  total  transmissivity  and  lateral  gradient, 
multiplied  by  the  width  of  flow,  W. 
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onn  non 


PROGRAM  "kin. f  " 


VARIABLES  AMD  PARAMETERS 


C 


* 


PARAMETER (DX=400  .  ,  DY=400.,  NROW=92,  NC0L=114,  M  =  3,  MW* 3, 
NCLASS=1 7 ,  iIN=10 ,  iRAIN=ll,  10RDER=12 ,  iHYDRO= 1 3 . 
iDIR=14,  iDISTS=15 ,  iDISTG=16,  iSLOPE*!',  iGAGE=13, 
iSC=l 9 ) 

REAL  CLASS (NRCW,  NCOL)  ,  ZWTm (NROW, NCOL)  ,  SLOPE (NRCW, MCOL) , 

ZF (NROW, MCOL) ,  ZT {NROW, MCOL) ,  NSTOm (NRCW, NCOL > , 

RQINm (NROW, NCOL)  ,  TIMEG (NROW, NCOL) ,  HYDRO (240), 

KOdat (NCLASS) ,  THSdat (NCLASS) ,  THRdat (NCLA33 ) ,  Edac ! MCLA33 ) 
Y  (N)  ,  YT (N)  ,  V ( NW )  ,  WK(87),  SS (NW) ,  LENGTH ( 3 )  ,  WIDTH (3), 

E,  F,  RR,  Ri,  RRF,  RAINRATE,  RQIN,  QOUT 
HR,  HI,  HRF,  AHR,  AHRF ,  NSTO,  MSTOold, 

T,  DT,  STEP,  RAINSTEP ,  TEND,  Told,  K,  S,  W,  WXY, 

ERROR,  HMIN,  TOL 

INTEGER  DISTS (NROW, NCOL) ,  DISTG (NROW, NCOL) ,  DIR (NRCW, MCCLi , 

GAGE (NRCW,  NCOL)  ,  PI  (8),  ?J(8),  IND,  IER,  ORDER,  £0 
CHARACTER* 8 0  DIRfile,  SLOPEfile,  DISTSfile,  DISTGf lie, 

SCfile,  CRDERfile,  GAGE file,  RAINfile,  HYDRO f lie, 

ZWTfile,  HEADER 

COMMON/REAL/RAIN  (20)  ,  nGAGE,  K0,  F,  E,  THS ,  THR,  P.i, 

NSTOold,  W,  S,  RQIN,  TO,  RAINRATE,  iTYPEl,  ZWT,  TRIGG 
EXTERNAL  MAINFCN,  ReFCN,  THFCN,  RAINFCN 
DATA  PI/0, 1, 1, 1, 0, -1, -1. -1/ 

DATA  PJ/-1, -1, 0, 1, 1, 1, 0, -1/ 

DATA  KOdat/  21.5  ,  3.6  ,  0.25  ,  0.25  ,  45. 

16.6  ,  21.8  ,  20.6  ,  20.6  ,  20.6  ,  20.6 


25 
21 , 


3.6 
.  3 


DATA  THSdat/  0.53  ,  0 

.52  , 

3.48  ,  3.48  , 

3.56 

,  3.56  ,  3.53  , 

z 

0.49  ,  0.52  ,  0 

.48  , 

3.52  ,  0.52  , 

3.52 

^  e  n 

z 

0.49  ,  0.48  ,  0 

.25  / 

DATA  THRdat/  0.02  ,  0 

.036  , 

3.09  ,  3.09 

9  ,  1  .  C  9  ,  *  .  1  C  9 

z 

0.109  ,  0.064  ,  0 

.036  , 

3.072  ,  3.07 

p  , 

3~2  ,  ;.3_  , 

i 

0.03  ,  0.041  ,  3. 

02  / 

DATA  Edac/  3.5  ,  3.6 

,  7.5 

,  7.5  ,  3.5  , 

-  C 

3.6,  2.6  ,  2  . 3  , 

3.4  ,  3.6  ,  3.6  , 

3.6  , 

3.6  ,  3.5  , 

J  .  C  t 

3  .  4 

READ 

INPUT 

5 

OPEN  <UNIT=iIN, F I LE  = ' k 

in .  m ' 

,  STATUS  =  ' OLD' 

) 

READ ( i IN, 10 ) SCfile,  DIRfile,  SLOPEfile,  DISTSfile, 
ORDERfile,  GAGEfile,  RAINfile,  ZWTfile,  HYDRO: 
FORMAT ( 9 (A8  0 / ) A8C ) 

READ  ( iIN,  *  )  F,  Ri,  Vo,  Vs,  iTYPE ,  .-.GAGE, 

RAINSTEP,  HMIN,  TOL,  TEND 

CLOSE ( i IN) 
iTYPEl  =  iTYPE 


OPEN (UNIT=iSC, FILE=SCf lie, 
OPEN  (L’NIT=iDIR,  FILE=D IRf  i  ' 
OPEN (ON IT  =  iD I STS, FILE =D  I  ST 
OPEN (ONIT=iDI3TG, FILE=DI3T 
OPEN (UNIT= iSLCPE , FILE=SLO? 
OPEN (UNIT=iZWT, FILE=ZWTf il 


STATUS =' 0L 
e, STATUS =' 
Siile, STAT 
jf 1*9,  S. AT 


Efile, STATUS-' 
e, STATUS =' OLD' 


READ (i SC, 20) HEADER 
READ (iDIR, 20) HEADER 
READ (iDISTS, 20) HEADER 
READ ( iDISTG, 20) HEADER 
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o  n  n 


READ  (iSLOPE,  20  HEADER 
READ (iZWT, 20) HEADER 
20  FORMAT (A80) 

DO  30  1  =  1 , NROW 

READ ( iSC, * )  (CLASS (I,  J)  ,  J=l,  NCOL) 
READ (iDIR, *)  (DIR ( I , J)  , J=l,  NCOL) 

READ (iDISTS,  *)  (DISTS (I, J)  , J=l, NCOL) 
READ (iDISTG,  *)  (DISTG ( I ,  J)  , J=1 , NCOL) 
READ ( iSLOPE,  *)  (SLOPE (I, J)  ,  J=1 , NCOL) 
READ ( iZWT, ')  ( ZWTrn (I,  J) , J=l, NCOL) 

30  CONTINUE 

CLOSE (iSC) 

CLOSE (iDIR) 

CLOSE (iDISTS) 

CLOSE (iDISTG) 

CLOSE (iSLOPE) 

CLOSE  (iZWT) 


PRELIMINARY  COMPUTATIONS 


WIDTH  OF  FLOW 


WXY  =  (DX*DY)  /  SQRT (DX**2+DY**2) 

WIDTH (I) =DY 

WIDTH (2) =WXY 

WIDTH (3) =DX 

WIDTH (4) =WXY 

WIDTH  (5) =DY 

WIDTH  (6) =WXY 

WIDTH (7) =DX 

WIDTH (8) =WXY 

DLXY  =  SQRT (DX**2*DY**2) 

LENGTH  (1) =DX 
LENGTH (2) =DLXY 
LENGTH (3) =DY 
LENGTH (4 ) =DLXY 
LENGTH  (5) =DX 
LENGTH (6) =DLXY 
LENGTH (7) =DY 
LENGTH (8) =DLXY 


C  IIME  OF  TRAVEL  TO  STREAMGAGE 


DO  40  1=1, NROW 
DO  40  J=1 , NCOL 

TIMES ( I, J)  =  float (DISTS (I, J) ) /Vo  -  float (DISTG ( I, J; )  Vs 


OPEN  INPUT  FILES 


IF ( i TYPE . EQ .2) THEN 

OPEN (UNIT=iGAGE, FILE=GAGEf ile, STATUS=' OLD' , READONLY: 
READ ( tGAGE,  * ) NX,  NY,  DDX ,  DDY 
DO  5C  1=1, NROW 

READ (iGAGE,  * )  (GAGE ( I ,  J) , U=l, NCOL) 

CONTINUE 
CLOSE (iGAGE) 
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OPEN (UNIT=iRAIN,  FILE=RAINf ile ,  STATUS  =  '  OLD' , READONLY) 
END  IF 


C - C 

C  TIME  LOOP  (9999)  C 

c - - - c 


* 


'» 


nSTEPS  =  INT (TEND/STEP) 

DO  9999  iT  =  1, nSTEPS 
TO  =  (i.T-1)  ‘STEP 

IF  (ABS ( INT (TO/RAINSTEP ) -TO/RAINSTEP )  . LT . IE-3 ) THEN 
CALL  RAINFCN ( iTYPE) 

END  IF 


c - C 

C  HILLSLOPE  PIXELS  LOOP  (8888)  C 

C - - - C 


OPEN (UNIT=iORDER,FILE=ORDERf ile,  STATUS  =  ' OLD' ) 
READ (iORDER, *) nHPIXELS , nSPIXELS 
nPIXELS  =  nHPIXELS +nSPIXELS 

DO  8888  IP=1, nHPIXELS 
READ (iORDER,  *) I,  J 

IF  (iTYPE.EQ.2)  RR  =  AMAX1  (RAIN  (GAGE  (I,  J)  )  ,  Ri) 


C - 

C  INITIALIZE  VARIABLES 
C - 


Y ( 1 )  =  ZF ( I , J) 

Y  (2 )  =  ZT  ( I ,  J) 

Y (3)  =  0. 

SS(1)  =  Y  ( 1 ) 

SS (2 )  =  Y (2) 

SS  (3)  =  Y  (3) 

SC  =  CLASS (I, J) 

KO  =  KOdat(SC) 

THS  =  THSdat(SC) 

THR=  THRdac ( SC) 

E  =  Edat (SC) 

S  =  SLOPE ( I , J) 

W  =  WIDTH (DIR (I, J) ) 

ZWT  =  ZWTm (I, J) 

NSTOold  =  NSTOm ( I , J) 

RQIN  =  RQINm ( I, J) 

T  =  TO 

DT  =  STEP 

alfa  =  ATAND(S) 

TRIGO  =  SIND (alfa) *COSD (alfa) 


C - 

C  UPDATE  VARIABLES 

C - 


IF (ZWT . EQ . 0 . ) THEN 
T  =  TO+STEP 
HR  =  RR*STEP 
HRF  =  HR 

AHR  =  AHR  +  HR/nPIXELS 
AHRF  =  AHRF  *  HR/nPIXELS 
VRF  =  HRF* (DX*DY) /1000. 

iHOUR  =  MAX  (NINT  (TO  •*■  ( STEP  /  2 )  +-T IMEG  ( I ,  J)  )  ,  1 ) 
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HYDRO (iHOUR)  =  HYDRO (iHOUR)  +  VRF  /  3600. 

GOTO  60 
END  IF 

70  CONTINUE 

CALL  DREBS  (MAINFCN,  Y,  T,  N,  6,  1 , 0 ,  DT , HMIN, TOL,  V,  SS, WK,  IER) 

IF (IER. NE.O) WRITE (*, *) 'ERROR' ,  IER, (V(IV) , IV=1,NW) 

IF (T . LT . (TO+STEP-HMIN) ) THEN 
DT  =  TO+STEP-T 
GOTO  70 
END  IF 

Y  (1)  =  AMIN1  (Y  (1)  ,  ZWT) 

Y (2)  =  AMIN1 (Y (2) , Y  (1)) 

Y(2)  =  AMAX1 (Y (2) ,0.) 

NSTO  =  AMIN1 (  (NSTOold+ (AMINl (RR,K0) +RQIN-Ri) 'STEP-Y  (3) )  , 
&  (  (THS-THR)  *Y  (1)  - (Ri/KO) *  * (1/E) * (THS-THR) * (E/F) * 

&  (EXP (F/E*Y  (1) )  -1. )  )  ) 

HI  =  NSTO  -  NSTOcld  -  Y ( 3 )  +  (Ri-RQIN) 'STEP 

HR  =  RR'STEP 

HRF  =  HR  -  HI 

RRF  =  HRF/STEP 

RI  =  AMAX1 (HI, 0 . ) /STEP 

VRF  =  HRF* (DX'DY) / 1000 . 

AHRF  =  AHRF  +  HRF/nPIXELS 
AHR  =  AHR  +  HR/nPIXELS 
AHI  =  AHI  +  HI/nPIXELS 

iHOUR  =  MAX (NINT(T0+ (STEP/2) +TIMEG(I, J) ), 1) 

HYDRO (iHOUR)  =  HYDRO (iHOUR)  -  VRF  /  3600. 


STORE  RESULTS 


ZF ( I , J)  =  Y(l) 

ZT ( I ,  J)  =  Y (2 ) 

NSTOm ( I , J)  =  NSTO 
RQOUT  =  Y (3) /STEP 

ID  =  DIR ( I , J) 

II  =  I  -  PI  (ID) 

JJ  =  J  -  P  J  ( ID) 

RQINm (II, JJ)  =  RQINm (II, JJ)  +  RQOUT 


60  CONTINUE 

RQINm ( I, J)  =  0. 


5888  CONTINUE 


3888 


STREAM  PIXELS  LOOP  (7777) 


00  7777  iSP  =  1 , nSP IXELS 
READ (iORDER, *) I, J 

IF ( iTYPE . EQ . 2 ) RR  =  PAIN (GAGE ( I , J) ) 
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c - 

C  INITIALIZE  VARIABLES 
C - 


RQIN  =  RQINm (I,  J) 
T  =  TO 


c - 

C  UPDATE  VARIABLES 
c - 


HR  =  RR*STEP 

RRF  =  RR  +  RQIN 

HRF  =  RRF  *  STEP 

VRF  =  HRF* (DX*DY) /1000. 

AHR  =  AHR  +  HR/nPIXELS 
AHRF  =  AHRF  +  HRF/nPIXELS 

iHOUR  =  MAX (NINT (T0+ (STEP/2) +TIMEG (I, J) ), 1) 
HYDRO (iHOUR)  =  HYDRO (iHOUR)  +  VRF  /  3600. 

RQINm (I, J)  =  0 


c - 7777 - C 

7777  CONTINUE 

CLOSE (iORDER) 


C - 9999 - C 

9999  CONTINUE 

CLOSE (iRAIN) 


c - 

C  WRITE  RESULTS 

c - 

OPEN (UNIT=iHYDRO,  FILE=HYDROf ile,  STATUS  =  ' NEW' ) 

WRITE (iHYDRO, *) 'END  TIME  [hr]  =  ',T 

WRITE (iHYDRO, *) 'CUM  RAIN  [mm] ='  , AHR, '  CUM  INFILT  [mmj=',AHI, 
&  '  CUM  RUNOFF  [ mm ] = ' , AHRF 

do  80  iT=l , TO+24 

80  WRITE ( iHYDRO, 90) iT*RAINSTEP,  HYDRO (iT) 

90  FORMAT (14 , F9 . 2) 

CLOSE (iHYDRO) 

END 


C 


& 

-  & 

& 

& 

Re  -  ReFCN (Y(1),Y(2),Y(3),T) 
THi  =  THFCN (Ri , Y ( 1 ) ) 

TH1  =  THFCN (Re, Y(l) ) 


SUBROUTINE  MAINFCN (N, T, Y, YT) 

PARAMETER (DX=400. ,  DY=400.) 

REAL  Y (N) ,  YT (N) 

REAL  K0,  F,  E,  THS ,  THR,  Ri,  RR,  RRF,  TOL,  HMIN, 

T,  STEP,  TEND,  HI,  NSTOoid,  AHR,  HU,  MS, 

TO,  AHRF,  HRF, 

W,  S  ,  V 

COMMON /REAL /RAIN  (20)  ,  nGAGE ,  K0,  F,  E,  THS,  THR,  Ri,  RR, 

NSTOoid,  W,  S,  RQIN,  TO,  RAINRATE,  iTYPEl,  ZWT,  TRIGG 


C- 

r* 

c- 


SUBROUTINE  "MAINFCN" 
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TH2  =  THFCN (Re, Y  (2)  ) 
Ks  =  K0*EXP (-F*Y  (1) ) 


SATURATION-EXCESS 


IF ( Y  { 1 )  .GE.ZWT) THEN 
YT (1)  =  0. 

YT ( 2 )  =  -  Re  /  (THS-TH2) 
ELSE 


UNSATURATED  FRONT 


IF (  (Y  (1)  -Y (2)  )  . LT . IE- 3 . AND . Ks .GT .  Re)  THEN 
YT ( 1 )  =  (Re-Ri)  /  (THl-THi) 

YT (2 )  =  YT (1 ) 

ELSE 


SATURATED  FRONT 


YT ( 1 )  =  (Ks-Ri)  /  (THS-THi) 

IF (THS-TH2 .GT . IE-3) THEN 
YT (2)  =  (  Ks  +  RLS  -  Re  )  /  (THS-TH2 ) 

ELSE 

IF (Ks .GT . Re) THEN 
SIGN  =  1. 

ELSE 

SIGN  =  -1. 

END  IF 

YT (2)  =  SIGN  *  E/ (THS-THR)  *  Re 
END  IF 

END  IF 

END  IF 

YT  ( 3 )  =  KO*TRIGO* (  (EXP(F*(Y(1)  -  AMAX1 (Y (2) , 0) > )  -  1.) 
&  /F-Y ( 1 ) +AMAX1 ( Y (2 )  ,  0 . )  ) 

&  /  ( EXP (F*Y{1)) *400000.) 

END 


FUNCTION  " ReFCN" 


REAL  FUNCTION  ReFCN (XI , X2 , X3 , X4 > 

C  XI : =Y ( 1 ) =ZF ,  X2=Y (2) =ZT,  X3=Y (3) =QOUT,  X4=T 

REAL  XI , X2 , X3, X4 , 

&  K0,  F,  E,  THS,  THR,  Ri,  RR, RQIN, NSTOold,  TO 

COMMON/REAL/RAIN (20)  ,  nGAGE ,  K0,  F,  E,  THS,  THR,  Ri,  RR, 

&  NSTOold,  W,  S,  RQIN,  TO,  RAINRATE,  iTYPEl,  ZWT,  TRIGO 

IF (iTYPEI .EQ. 1) THEN 

ReFCN  =  RR 

ELSE 


IF ( (F/E*X2) .LT. IE-4) THEN 

ReFCN  =  RR  !  good  assumption  if  ZT  is  very  small! 
ELSE 


& 

& 

& 


ReFCN  =  K0M  (  NSTOold  +  (RR+RQIN)  *  (X4-T0 )  -  X3 

+  (Ri/KO) ** (1/E) * (THS-THR)  *  (E/F) * (EXP (F/E*X1) -1  .  ) 
-  (THS-THR) » (X1-X2)  )  / 

(  (THS-THR) * (E/F) * (EXP (F/E*X2) -1 . )  )  )**E 


END  IF 
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n  n  n 


mood  t-'ono  non  n  nnn 


END IF  * 

END 

r 


FUNCTION  "THFCN" 


REAL  FUNCTION  THFCN (XI, X2) 

Xl=Re,  X2=ZF  or  ZT 
REAL  XI, X2, 

&  KO,  F,  E, THS, THR, Ri, RR, RQIN, NSTOold, TO 

COMMON/REAL/RAIN (20) ,  nGAGE,  KO,  F,  E,  THS,  THR,  Ri,  RR, 

&  NSTOold,  W,  S,  RQIN,  TO,  RAINRATE,  iTYPEl ,  ZWT,  TRIGO 

THFCN  =  AMIN1 (  (  (X1/K0) ** (1/E) *EXP (F/E*X2) * (THS-THR) +THR  ), 
&  THS  ) 


END 


SUBROUTINE  "RAINFCN” 


SUBROUTINE  RAINFCN ( iTYPE) 

COMMON/REAL/RAIN (20)  ,  nGAGE,  KO,  F,  E,  THS,  THR,  Ri ,  RR, 

&  NSTOold,  W,  S,  RQIN,  TO,  RAINRATE,  iTYPEl,  ZWT,  TRIGO 

GOTO (1,2) iTYPE 


UNIFORM  RAINFALL 


CONTINUE 
RR  =  RAINRATE 
RETURN 


OBSERVED  RAINFALL 


CONTINUE 

READ (11,*)  (RAIN (iG) , iG=l,  nGAGE) 
RETURN 

END 


\ 
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